2.1 Was ist R?

R ist eine Programmiersprache, die zur Datenanalyse entwickelt wurde und damit eine Alternative zu Programmen wie SPSS, SAS und den Statistik-Funktionen von Excel. Die zwei größten Vorteile von R gegenüber point-and-click Programmen sind:

1. Der Funktionsumfang ist grundsätzlich nicht beschränkt → durch Programmierung unendlich erweiterbar → Über 10.000 Software-Pakete, die verschiedenste Analyseverfahren ermöglichen
2. Die Analysen mit R sind reproduzierbar und transparent! Analyse kann ohne Aufwand wiederholt werden, z.B. nach Fehler in den Rohdaten. Jeder, der den Code sieht (und versteht), kann nachvollziehen, wie die Daten verwendet worden sind.

Weitere Vorteile sind:

• die Software ist open source also kostenfrei verfügbar
• sie läuft auf allen gängigen Betriebssystemen
• es gibt verschiedene Editoren, die die Nutzung erleichtern, der bekannteste ist R-Studio
• es gibt eine sehr aktive Community im Netz, die Hilfe bietet, oft innerhalb von Minuten bis Stunden

R ist auch eine Investition in Ihre (berufliche) Zukunft. Außer Daten auswerten kann man damit

• komplexe und schöne Graphiken erzeugen (siehe Abbildung unten)
• Berichte, Paper und Bücher schreiben (auch Bachelor- und Masterarbeiten)
• Web-Seiten erstellen (wie diese hier)
• Apps entwickeln
• Simulationsmodelle implementieren
• Eigene Methoden entwickeln und veröffentlichen

Beispiele für R-Graphiken

Ganz allgemein gilt R zunehmend als die Standardsprache für statistische Problemstellungen sowohl in der Wirtschaft als auch in der Wissenschaft (Amirtha et al., 2014).

2.2 R und R-Studio installieren

1. Zuerst sollten Sie die eigentliche Programmiersprache R installieren. Auf dieser Seite finden Sie Downloads für verschiedene Betriebssysteme:

https://www.r-project.org

Sie gelangen hier auf eine Seite mit einem Link zum R-download. Folgen Sie diesem Link und suchen Sie dann einen CRAN-mirror aus. Im Prinzip ist egal, welchen Sie wählen (auf allen sind die gleichen Ressources verfügbar, deshalb ‘Spiegel’), es macht aber Sinn, einen Mirror in der Nähe zu wählen. Es wird dann ein Installationsprogramm (z.B. eine exe-Datei oder eine .deb Datei) heruntergeladen. Wenn Sie diese öffnen, werden Sie durch die Installation geführt - bei den Auswahlmöglichkeiten passen normalerweise voreingestellten Optionen.

2. Dann (und erst dann! sonst kann es zu Konfigurationsproblemen kommen) installieren Sie auch noch einen Editor für R. Der Editor ist hauptsächlich dafür da, den Code, zum Beispiel für statistische Analysen, zu schreiben, zu kommentieren und zu speichern. R-Studio ist der bekannteste und wahrscheinlich beste Editor:

https://www.rstudio.com/products/rstudio/download/#download

Wenn Sie R-Studio öffnen, sehen Sie vier Fenster: den eigentlich Editor, in dem Sie R-code schreiben, die R Konsole, in der die aktive R-Session läuft, eine Übersicht über den Arbeitsspeicher und unten rechts ein Fenster in dem je nach Aktion Plots, Hilfedateien oder die Ordnerstruktur und Dateien des Verzeichnisses, in dem der R-Prozess gestartet wurde, zu sehen sind. Eventuell fehlt auch der Editor oben links noch, das ändert sich aber, wenn Sie ein neues Skript anlegen (siehe ‘Erste Schritte mit R’)

Ansicht von R-Studio mit 4 Fenstern

2.3 Erste Schritte mit R

• Ein neues Skript im Editor öffnen, um Code schreiben zu können (File → new file → R-Script)
• R als Taschenrechner benutzen: In den Editor

5 + 7

schreiben → auf den ‘Run’ Button klicken (über dem Editor) → der Code wird an den laufenden R Prozess in der Konsole geschickt und das Ergebnis in der Konsole ausgegeben

• Eine Variable anlegen: im Editor

Ergebnis1 <- 5 + 7

eingeben → Run-Button anklicken → nichts passiert außer dass die Zeile in der Konsole auftaucht? → das Ergebnis wird in der Variable Ergebnis1 gespeichert, diese Variable wird jetzt auch in der Übersicht des Arbeitsspeichers gelistet. Wenn man in der Konsole

Ergebnis1

eingibt, wird der Wert der Variable ausgegeben. Anstelle des Pfeils kann auch ein = verwendet werden.

• Kommentare schreiben: Damit andere oder unsere zukünftigen Ichs später noch wissen, was wir warum gemacht haben, kann und sollte man in das R-Skript Kommentare schreiben. Wenn solche Zeilen zum R-Prozess geschickt werden, werden sie einfach ignoriert. Ein Kommentar wird mit einem Raute-Zeichen (hash) eingeleitet.

# Hier benutzen wir R als Taschenrechner
Ergebnis1 <- 5 + 7

• das Skript speichern: File → Save Dabei werden nicht die Variablen, Ergebnisse und Daten gespeichert (das lernen wir später). Sie können aber mit einem Knopfdruck das gesamte Skript laufen lassen und so alle Arbeitsschritte wiederholen → alle Ergebnisse wieder herstellen.

2.4 Importieren von Daten aus Excel/LibreOffice nach R

1. Daten in Excel/LibreOffice eingeben
   • Häufige Fehlerquelle! Am besten zu zweit: einer liest vor, einer tippt. Zwei/drei durch zwo/drei disambiguieren. Wenn nur eine Person die Daten eingibt, am besten den Nummernblock der Tastatur verwenden und Zahlen blind eingeben, also nur auf Zettel und Bildschirm schauen. Eingegebene Daten danach nochmal prüfen.
   • Keine leeren Zellen im Datenblock → wenn die Werte fehlen: “NA” (= not available)
   • Erste Reihe wird per default als Spaltenüberschrift interpretiert, hier aussagekräftige aber nicht zu lange Überschriften wählen
   • Zeichenketten (wie Überschriften) sollten aus einem einzigen Wort bestehen. Wenn nötig, zum Beispiel Unterstriche verwenden um einzelne Wörter zu verbinden
   • Als Dezimaltrenner sollte ein Punkt genutzt werden → “7.5” nicht “7,5”
     
   • Pro Beobachtung sollte es eine Zeile geben
     
2. Daten als csv Datei speichern: ‘Save as’ -> Dateityp aussuchen (csv, Text csv oder ähnlich) -> Character set sollte UTF-8 sein und der Field delimiter (Trennzeichen) ein Komma
3. RStudio öffenen und ein neues Projekt anlegen: ‘File’ (ganz oben links) -> ‘New Project’ -> ‘Existing Directory’ -> ‘Browse’ -> den Ordner aussuchen, indem die Kirsch-Datei gespeichert ist -> ‘Choose’ -> ‘Create Project’
4. Ein neues Skript anlegen: ‘File’ -> ‘New file’ -> ‘R Script’
5. In diesem neuen Skript kirschen <- read.csv("Kirschen.csv") eingeben um die Daten zu importieren.

Damit Sie erstmal keine Daten eingeben müssen, können Sie die Kirsch-Daten hier herunterladen und dann mit Punkt 3 weitermachen:
https://ecampus.uni-bonn.de/goto_ecampus_file_3525090_download.html

2.5 Daten kontrollieren

Wenn Sie die Daten in R eingelesen haben, sollten Sie als erstes kontrollieren, ob der Import geklappt hat. Eine erste Übersicht bekommen Sie mit den Funktionen summary() und head(). Der Datensatz heißt ‘kirschen’, weil wir das beim Import so angegeben haben.

summary(kirschen)
# gibt eine Zusammenfassung der Daten mit Mittelwerten etc. aus

head(kirschen)
# Zeigt die ersten 6 Zeilen des Datensatzes, wie er in R eingelesen wurde

Hier prüft man unter anderem, ob die Daten den richtigen Typ haben, also Zahlen mit denen wir rechnen wollen ‘numeric’ sind (dann werden ein paar Statistiken, wie Minimum, Mittelwert etc angegeben) und Kürzel für Versuchsbehandlungen Faktoren (dann sieht man, wie viele Beobachtungen es pro Gruppe = Level gibt). Ist das nicht der Fall, kann es jetzt geändert werden. Hier sind die Untersuchungsgruppen zum Beispiel als ‘character’ eingelesen worden (reine Zeichenketten)deshalb wandeln wir sie mit der Funktion ‘as.factor’ in Faktoren um:

kirschen$Unterlagenklon <- as.factor(kirschen$Unterlagenklon)
kirschen$Inokulat <- as.factor(kirschen$Inokulat)

Auf die einzelnen Variablen können Sie zugreifen, indem Sie zuerst den Namen des Datensatzes (= ‘kirschen’) verwenden und dann nach dem Dollar Zeichen die Spaltenüberschrift der Variable.

Auch einfache Plots sind ein gutes Mittel, um eine erste Vorstellung der Daten zu bekommen und Auffälligkeiten wie extreme Werte (verrutscher Punkt in den Zahlen) zu entdecken. Dafür eignen sich

• die einfach plot-Funktion

  plot(kirschen$Inokulat, kirschen$Trockengewicht, 
   main = "Inokulat", ylab = "Trockengewicht der Blätter [g]")

• oder eine schönere Version mit dem Paket ‘ggplot2’. Dazu muss das Paket ‘ggplot2’ erst installiert werden. In R-Studio geht das über ‘Tools -> Install Packages’ -> a search in the ‘Packages’ field.

  library(ggplot2)    
  ggplot(kirschen, aes(x=Inokulat, y=Trockengewicht)) +  
    geom_boxplot(notch=TRUE)

2.6 Hilfe bekommen

Hilfe zu finden ist zu Beginn des Arbeitens mit R die wichtigste Kompetenz überhaupt. Am einfachsten ist es natürlich, wenn Sie direkte Unterstützung durch jemanden bekommen, der sich mit R auskennt (und das ist bei der Bachelor-Arbeit - zumindest im Gartenbau - natürlich der Fall). Es gibt aber auch eine Reihe anderer Möglichkeiten:

1. in R selbst
   • ein “?” vor den Funktionsnamen setzen ?sd() → es öffnet sich eine Hilfedatei,
   • help.start() → öffnet Hilfeseite mit wichtigen Dokumenten
2. im Internet
   • Bei der Suche in der Suchmaschine ein R zu der Frage oder der Fehlermeldung hinzufügen (Fehlermeldungen gründlich lesen!)
   • Nach Tutorials und online Büchern suchen, zum Beispiel YaRrr - the pirates guide to R https://bookdown.org/ndphillips/YaRrr/
   • Cheat cheets
     …bieten einen Überblick der verfügbaren Funktionen und kurze Beschreibungen zu einem bestimmten Thema. Hier finden Sie eine Sammlung von cheat cheets: https://www.rstudio.com/resources/cheatsheets/ Hier ein cheat sheet mit den grundlegenden Funktionen: https://iqss.github.io/dss-workshops/R/Rintro/base-r-cheat-sheet.pdf

3.1 R-Projekte

R-Projekte sind gewissermaßen das zu Hause von allen R-Skripten und Daten, die zu einem Projekt (zum Beispiel einem Experiment) gehören. Dabei ist ein Projekt aber mehr als nur ein Ordner. Die Vorteile eines Projekts im Vergleich zu einer einzelnen neuen R Script-Datei sind, dass

• beim Öffnen des Projekts das Arbeitsverzeichnis automatisch auf den Ordner gesetzt wird, in dem sich die Projekt-Datei befindet. Man kann dann Dateien direkt aufrufen/einlesen ohne den vollständigen Pfad angeben zu müssen (höchstens Unterordner, wenn man denn welche angelegt hat). In Bezug auf das letzte Kapitel sollte jetzt also kirschen <- read.csv("Kirschen.csv") immer funktionieren, wenn die Datei ‘Kirschen.csv’ im gleicher Ordner, wie die .Rproj-Datei liegt.
• RStudio bei einem Neustart alle Dateien wieder öffnet, die beim letzten Schließen offen waren
• man gleichzeitig mehrere Projekte offen haben kann, ohne durcheinander zu kommen
• man das ganze Projekt komplett mit richtigen Verknüpfungen zwischen Daten/Dateien weiterreichen kann und die Skripte (zumindest bei gleicher R-Version) auf jedem Rechner reproduziert werden können.

Sie können ein RStudio-Projekt sowohl in einem neuen Verzeichnis als auch in einem vorhandenen Verzeichnis, in dem Sie bereits R-Code und Daten haben, erstellen. Um ein neues Projekt zu erstellen, verwenden Sie den Befehl ‘File -> New Project’.

Es öffnet sich dann ein Fenster, in dem Sie entscheiden können, ob Sie das Projekt in einem neuen Verzeichnis, einem schon existierenden Verzeichnis oder mit Versionskontrolle (zum Beispiel mit git) anlegen wollen. Letzteres ist sehr sinnvoll, würde hier aber den Rahmen sprengen.

Durch das Anlegen eines neuen Projekts in RStudio wird eine Projektdatei (mit der Erweiterung .Rproj) im Projektverzeichnis erstellt. Diese Datei kann später als Verknüpfung zum direkten Öffnen des Projekts verwendet werden.

3.2 Datentypen und -transformationen

Das Einlesen der Daten sollte jetzt kein Problem mehr sein, allerdings haben die Daten manchmal nicht den Datentyp der für bestimmte statistische Auswertungen oder Abbildungen gebraucht wird. Zuerst eine Übersicht über die gängigsten Typen:

Mit der Funktion str()kann man testen, als welche Datentypen R die Daten interpretiert und bekommt eine Übersicht der Struktur des Datensatzes:

str(kirschen)

Wie oben erwähnt, werden für einige Auswertungen/Abbildungen manchmal andere Datentypen benötigt, als R erkannt hat. So werden zum Beispiel manchmal die Bezeichnungen für Behandlungsgruppen als ‘character’ eingelesen, für eine Varianzanalyse werden aber Faktoren benötigt. In diesem Fall kann der Datentyp mit der Funktion as.factor() umgewandelt werden.

kirschen$Inokulat <- as.factor(kirschen$Inokulat)

Das geht natürlich nur, wenn der neue Datentyp zu der umzuwandelnden Variablen passt: ‘Gi_12’ kann nicht in eine numerische Variable umgewandelt werden, wohl aber die Zahlen 1, 2, 3, .. in Faktoren (die dann als “1”, “2”, “3” ausgegeben werden). Andere häufig genutzte Funktionen sind as.numeric() , as.character() und as.Date().

3.3 Umgang mit NAs

Auch NAs können manchmal zu Fehlermeldungen oder unvollständigen Abbildungen führen. Da NA nicht Null bedeutet, sondern eher ‘wissen wir nicht’, kann zum Beispiel nicht ohne Weiteres ein Mittelwert berechnet werden.

Daten <- c(5, 3, 9, NA, 4, NA)
mean(Daten)

Die meisten Funktionen haben aber die Option, NA-Werte bei der Berechnung auszuschließen:

mean(Daten, na.rm = TRUE)

Schauen Sie also am besten zuerst in der Hilfe zu der Funktion , die Sie verwenden möchten (?mean), ob es ein Funktionsargument gibt, das den Umgang mit NA-Werten festlegt und ändern Sie die Einstellung entsprechend. Wenn es gar nicht anders geht, kann man als Alternative auch die Zeilen in den NAs vorkommen vor der Analyse herausfiltern. Wie solche Daten-Operationen funktionieren, lernen Sie im folgenden Kapitel.

3.4 Datensätze umstrukturieren

Manchmal ist es nötig, den eingelesenen Datensatz umzustrukturieren, Teildatensätze herauszufiltern oder nue Variablen zu berechnen. Dazu sollte man am besten das Paketbündel ‘tidyverse’ installieren. Hat man das getan, steht auch eine weitere Konvertierungsfunktion zur Verfügung, die aus einem normalen Datensatz (data.frame) ein ‘tibble’ macht. Tibble haben ein paar Vorteile, zum Beispiel, dass standardmäßig nur die ersten 10 Zeilen ausgegeben werden und zu jeder Spalte der Datentyp gleich dabei steht (wobei Fließkommazahlen hier nicht als numeric sondern double (dbl) bezeichnet werden. ‘fct’ steht für Faktor, ‘int’ für Integer).

kirschen <- as_tibble(kirschen)
kirschen

Daten <- c(5, 9, 2, 5, 3, 9, 5, 3, 5)
Ergebnis <- round(sqrt(sum(Daten^2)), 3)
Ergebnis

Daten_quad <- Daten^2
Daten_quad_sum <- sum(Daten_quad)
Daten_quad_sum_wurzel <- sqrt(Daten_quad_sum)
Ergebnis <- round(Daten_quad_sum_wurzel, 3)

Ergebnis <- Daten^2 %>%   
            sum() %>%     
            sqrt() %>%    
            round(3)      

library(tidyverse)
# mit der Funktion 'tibble()' legt man einen Datensatz an
Daten1 <- tibble(spalte1 = c(1, NA, 5, 7), spalte2 = c("a", "b", NA, "d"))
# wir schauen uns den Datensatz an:
Daten1
# wir löschen die Zeilen mit NA heraus
Daten2 <- Daten1 %>% drop_na()
# uns sehen uns das Ergebnis an:
Daten2

kirschen %>% select(Trockengewicht, Blattfläche) %>% head()

kirschen %>% select(-Trockengewicht, -Blattfläche) %>% head()

kirschen_neu <- kirschen %>% select(-Trockengewicht, -Blattfläche)
head(kirschen_neu)

kirschen_nurPDV <- kirschen %>% filter(Inokulat == "PDV")
summary(kirschen_nurPDV)

plot(kirschen_nurPDV$Inokulat, kirschen_nurPDV$Frischgewicht)

kirschen_nurPDV <- kirschen %>% 
                   filter(Inokulat == "PDV") %>%
                   droplevels()
                   
 plot(kirschen_nurPDV$Inokulat, kirschen_nurPDV$Frischgewicht)                  

kirschen %>% mutate(Gesamttrieblänge_m = Gesamttrieblänge/100)

Wachstum <- tibble(Baum = c("Nr1", "Nr2", "Nr3", "Nr4"), 
                   Zeitp1 = c(338,459,558,439), 
                   Zeitp2 = c(437,560,729,658), 
                   Zeitp3 = c(469,579,937,774))
                   
Wachstum_lang <- Wachstum %>% 
                 pivot_longer(
                 cols = c(Zeitp1, Zeitp2, Zeitp3),
                 names_to = "Zeitpunkt",
                 values_to = "Größe_cm")
                 
# cols: die Zeilen die zusammengefügt werden sollen
# names_to: neuer Spaltenname des Faktors
# values_to: neuer Spaltenname der Daten

# 'Zeitpunkt' und 'Baum' wurden noch nicht als Faktor erkannt, deshalb:
Wachstum_lang$Zeitpunkt <- as.factor(Wachstum_lang$Zeitpunkt)
Wachstum_lang$Baum <- as.factor(Wachstum_lang$Baum)

ggplot(Wachstum_lang, aes(x = Zeitpunkt, y = Größe_cm, color=Baum)) + 
   geom_point() +   labs(y="Größe (cm)", x="Zeitpunkt")

Wachstum_ww <- Wachstum_lang %>% 
               pivot_wider(
               names_from = "Zeitpunkt",
               values_from = "Größe_cm"
               )

3.5 Kontinuierliche Daten auswerten

Bevor wir in die eigentliche Datenanalyse einsteigen, ist es wichtig, dass Sie das Grundprinzip der statistischen Tests verstanden haben. Auf die Gefahr hin, dass es eine Wiederholung für einige von Ihnen ist, wollen wir deshalb nochmal die Idee der Varianzanalyse, die den meisten Tests zugrunde liegt, erklären. Am Ende dieses Kapitels

• haben Sie die Logik der Varianzanalyse verstanden
• können Sie etwas mit den Begriffen ‘Occams Razor’ und ‘Homoskedastizität’ anfangen
• sind Sie in der Lage, eine einfache Varianzanalyse mit kontinuierlichen Daten in R durchführen
• wissen Sie, wie man die Voraussetzungen einer ANOVA visuell testet
• wissen Sie auch, dass nicht die Daten als Gesamtheit normalverteilt sein müssen, sondern die Daten innerhalb der Behandlungsgruppen, bzw. die Residuen (Fehler, Abstände zum Mittelwert der Gruppe).

3.6 Das Grundprinzip der Varianzanalyse

Wir nehmen an, wir haben einen sehr einfachen Versuch mit Tomaten durchgeführt, in dem die Hälfte aller 40 Pflanzen gedüngt wird und die andere Hälfte als Kontrollgruppe ungedüngt bleibt. Am Ende bestimmen wir das Gesamtgewicht der Tomaten pro Pflanze als ein Maß für den Ertrag. Jetzt interessiert uns, ob die Düngung einen signifikanten Effekt auf den Ertrag hat. Unsere Hypothese ist: Die Düngung hat einen Effekt auf den Ertrag. Statistisch ausgedrückt bedeutet das, dass die Ertragsdaten aus zwei Normalverteilungen mit unterschiedlichen Mittelwerten stammen. Diese Hypothese ist auf folgender Abbildung auf der linken Seite dargestellt. Das zugehörige statistische Modell hat zwei Parameter: a ist der Mittelwert der Normalverteilung ohne Dünger, b die Differenz zwischen a und dem Mittelwert für die Normalverteilung mit Dünger (die Effektgröße). Die Variable Dünger nimmt entweder den Wert 0 (ohne Dünger) oder 1 (mit Dünger) an. Die Streuung der Daten (sie liegen nicht alle genau auf dem Mittelwert) wird durch den Fehlerterm aufgefangen.

Hypothesen, statistische Modelle und Umsetzung in R. Links: Düngung hat einen Effekt auf den Ertrag -> die Daten stammen aus zwei Normalverteilungen mit unterschiedlichen Mittelwerten; rechts (Nullhypothese): Düngung hat keinen Effekt auf den Ertrag -> die Daten stammen aus einer einzigen Normalverteilung.

Die Null-Hypothese zu diesem Modell ist: Die Düngung hat keinen Effekt auf den Ertrag -> die Daten stammen alle aus einer einzigen Normalverteilung mit dem Mittelwert a (rechte Seite der Abbildung). Entsprechend einfacher ist das statistische Modell.

In R formulieren wir diese beiden konkurrierenden Modelle mit der Funktion lm(). lm steht für ‘linear model’ und es schätzt Mittelwert (oder in anderen Fällen auch Achsenabschnitt und Steigung einer Regressionsgeraden), indem es die Werte so wählt, dass die Fehlerquadrate (quadrierter Abstand der Messungen zu den Mittelwerten) möglichst klein sind.

# hier lesen wir keine Daten aus Excel ein sondern geben sie einfach händisch direkt 
# in R ein. Zuerst Ertragsdaten von 10 ungedüngten Pflanzen
ohne <- c(417.1,558.8,519.1,611.7,450.3,461.9,517.3,453.2,533.1,514.7)
# dann von 10 gedüngten Pflanzen
mit <- c(581.9,517.4,641.3,659.8,587.8,483.3,703.0,589.5,532.4,569.3)
# hier werden die beiden Datensätze mit c() wie combine
#  zu einem Vektor (= Spalte) zusammengeführt 
Ertrag <- c(ohne, mit)
# jetzt generieren wir noch einen Vektor der die Behandlung anzeigt
Düngung <- gl(2, 10, 20, labels = c("ohne","mit"))
# Mit diesen beiden Vektoren können wir das statistische Modell formulieren, 
# dass der Ertrag von der Düngung abhängt. In R wird das mit der Tilde ~ ausgedrückt
#  (oben links auf der Tastatur). 
# Wir nennen dieses Modell 'model1'
model1 <- lm(Ertrag ~ Düngung)
# mit der folgenden Zeile schauen wir uns eine summary (Zusammenfassung)
#  des gefitteten Modells an

summary(model1)

Interpretation der Ausgabe einer model summary

In der Zeile ‘Intercept’ finden wir unter ‘Estimate’ den geschätzten Mittelwert für die Gruppe ohne Düngung, also das a in der Abbildung oben, linke Seite. Darunter, in der Zeile ‘Düngungmit’ ist die Differenz (das b aus der Abbildung) das auf das a aufaddiert werden muss, um die Schätzung für den Mittelwert der gedüngten Pflanzen zu erhalten.

Als nächstes formulieren wir die Nullhypothese

null_model <- lm(Ertrag ~ 1)
summary(null_model)

## 
## Call:
## lm(formula = Ertrag ~ 1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -128.04  -38.30  -12.39   43.08  157.85 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   545.15      16.68   32.69   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 74.59 on 19 degrees of freedom

Wie erwartet, erhalten wir jetzt nur einen Mittelwert über alle Tomatenpflanzen.

Welches der beiden Modelle ist nun besser? Das Modell mit Düngung als erklärender Variable passt sich besser an die Daten an -> kleinere Summe der Fehlerquadrate, ist aber komplizierter.

3.7 Das Prinzip der Parsimonie

Die Grundlagen für das Prinzip der Parsimonie, auch bezeichnet als Occam’s Rasor. Verändert nach einer Folie von Frank Schurr

Um zu entscheiden, ob eine erklärende Variable einen signifikanten Effekt hat, stellen wir die Frage: Wie wahrscheinlich ist es, die Daten so - oder mit noch größeren Unterschieden zwischen den Behandlungen - zu beobachten, wenn in Wahrheit die erklärende Variable keinen Effekt hat (= die Nullhypothese wahr ist bzw. alle Daten aus einer einzigen Normalverteilung stammen). Um diese Wahrscheinlichkeit ausrechnen zu können, müssen die Daten innerhalb der Behandlungsgruppen normalverteilt sein und die Streuung (Varianz) muss in etwa gleich sein. Das sind die wichtigsten Vorraussetzungen für eine Varianzanalyse, weil sonst evtl. die berechnete Wahrscheinlichkeit nicht stimmt. In R nutzten wir die Funktion anova() um die beiden Modelle zu vergleichen und die oben genannte Wahrscheinlichkeit auszurechnen:

anova(model1, null_model)

Interpretation der Ausgabe einer ANOVA

Die Wahrscheinlichkeit nach der wir suchen, finden wir unter Pr(>F) und ist der berühmte p-Wert. Die Wahrscheinlichkeit zufällig solche (oder noch extremere) Daten zu finden, wie wir sie für die Tomaten gefunden haben, obwohl die Düngung in Wahrheit keinen Einfluss auf den Ertrag hat, ist nur 0,86% und damit so unwahrscheinlich, dass wir den Effekt der Düngung als eindeutig signifikant interpretieren können. Damit ist model1 das bessere Model. Eine allgemeine, wenn auch etwas arbiträre Übereinkunft, ist, dass alle Variablen mit einem p-Wert von unter 5% bzw p < 0.05 als signifikant gelten.

Kritik an der Verwendung des p-Wertes
In den letzten Jahren gibt es mehr und mehr Kritik an dem starken Fokus auf den p-Wert. Zum einen, weil zumeist einfach das Signifikanzniveau verwendet wird, was fast immer verwendet wird, ohne weiter darüber nachzudenken: 5%. Dieses Niveau stellt einen Kompromiss zwischen der Fähigkeit eine Entdeckung zu machen und unserer Bereitschaft ein fehlerhaftes Testergebnis zu akzeptieren dar. Es sollte eigentlich einleuchtend sein, dass das Signifikanzniveau für verschiedene Fragestellungen auch unterschiedlich angesetzt werden sollte (was der Begründer der frequentistischen Statistik, Ronald Fisher auch schon 1956 selbst so geraten hat). Zum anderen ist der Nachweis eines signifikanten Effekts außer von der Effektgröße stark vom Informationsgehalt der Daten abhängig (wie groß ist die Stichproble). So kann es sein, dass, wenn man hunderte von Proben hat, eine Variable als signifikant nachgewiesen werden kann, ihr Effekt und die wissenschaftliche Relevanz aber nur minimal ist. Es ist also immer, immer wichtig, den Bezug zur eigentlichen Frage und zum untersuchten System zu behalten und die Ergebnisse mit gesundem Menschen- und Sachverstand zu interpretieren.

3.8 Annahmen der ANOVA überprüfen

Wie oben schon kurz angesprochen, müssen für glaubwürdige Ergebnisse einer Varianzanalyse einige Vorraussetzungen erfüllt sein

1. Unabhängigkeit der Messungen - das sollte beim Design des Experiments bereits bedacht worden sein.
2. Keine großen Ausreißer in den Daten - das geht am einfachsten durch das Plotten der Daten und evtl. Ausreißer herausnehmen.
3. Die abhängige Variable ist innerhalb der Behandlungsgruppen normalverteilt.
4. Die Varianz in jeder Gruppe sollte annähernd gleich sein. Das wird auch mit dem Begriff Homoskedastizität beschrieben.

Um die Normalverteilung der Daten innerhalb der Behandlungsgruppen zu prüfen, ist es einfacher, sich die Daten nicht gruppenweise anzuschauen, sondern einfach die Residuen/Fehler aller Daten (also den Abstand jedes Datenpunktes zum geschätzten Gruppen-Mittelwert) zu nehmen. Es gibt numerische Tests, um die Normalverteilung der Residuen zu prüfen, in den letzten Jahren setzt sich aber mehr und mehr eine visuelle Überprüfung durch. Am besten eignet sich dafür ein qq-Plot (Quantile-Quantile-Plot). Wenn die Punkte annähernd auf einer Geraden liegen, sind die Residuen normalverteilt. Da ‘annähernd’ ein dehnbarer Begriff ist, liefert die Funktion qqPlot aus dem Package ‘car’ ein Konfidenzintervall. Liegen die Punkte innerhalb dieses Intervalls, kann man von einer Normalverteilung ausgehen. Das Paket ‘car’ müssen Sie in R-Studio vorher installieren (Tools -> Install Packages) Im jupiter-Notebook sollten Sie zu Beginn der Zelle diesen Code einfügen

# Schnelles Repository: Posit Public Package Manager
options(repos = c(PPM = "https://packagemanager.posit.co/cran/__linux__/jammy/latest"))

# Hier die benötigten Pakete auflisten
pkgs <- c("car")

# Nur fehlende Pakete installieren
to_install <- setdiff(pkgs, rownames(installed.packages()))
if (length(to_install)) install.packages(to_install)

Dann geht es weiter:

library(car)

## Loading required package: carData

## 
## Attaching package: 'car'

## The following object is masked from 'package:dplyr':
## 
##     recode

## The following object is masked from 'package:purrr':
## 
##     some

qqPlot(resid(model1))

## [1] 17  4

Außer dem Plot gibt die Funktion die ‘Namen’ der Punkte aus, die im Bezug auf die x-Achse outlier darstellen (also Punkte, die besonders stark abweidchen,hier Punkt 4 und Punkt 17). Es lohnt sich, diese Punkte genauer anzuschauen (waren die Versuchsbedingungen hier aus irgendeinem Grund anders?) und evtl. auszuschließen, insbesondere, wenn sie außerdem außerhalb der Konfidenzlinien liegen.

Eine gute Erklärung wie qq-plots Funktionieren finden Sie hier: https://www.youtube.com/watch?v=okjYjClSjOg … auf Englisch und mit toller Intro :)

Um die Homoskedastizität visuell zu überprüfen, plotten Sie die Residuen gegen die von der Funktion lm() geschätzten Mittelwerte. Wenn die Daten der Gruppen etwa einen gleichen Bereich auf der y-Achse umfassen, sind die Varianzen homogen. Hierfür gibt es leider derzeit keine Konfidenzintervalle, eine Daumenregel ist aber, dass die Varianzen als homogen gelten, wenn die größte Varianz nicht mehr als 1,5fach größer ist als die kleinste Varianz. In unserem Beispiel reichen die Residuen beider Gruppen von etwa -100 bis 100, die Varianzen sind also annähernd gleich groß. Formal kann Varianzhomogenität mit dem Leven’s Test überprüft werden. R bietet dazu die Funktion leveneTest() aus dem package car. Wer interessiert daran ist, kann sich mit ?leveneTest die Hilfeseite mit Erklärungen zur Anwendung der Funktion anschauen.

plot(fitted(model1), resid(model1))
abline(h=0, lty = 3, col = "red")

Für unsere Analyse sind alle Annahmen erfüllt. Wir können also davon ausgehen, dass die Düngung in unserem hypothetischen Experiment tatsächlich einen signifikanten Effekt auf den Ertrag hat. In einer Arbeit würden wir den F-Wert (aus summary(model1)), die Anzahl der Freiheitsgrade und den p-Wert angeben, wobei p-Werte, die kleiner als 0.05 sind, üblicherweise mit p < 0.05 angegeben werden.

4.1 Nicht-normalverteilte Daten

Im letzten Kapitel haben wir gesehen, dass eine Voraussetzung für die Durchführung einer ANOVA die Normalverteilung der Residuen (Fehler, Abstände der Messpunkte zum Mittelpunkt) ist. Bei kontinuierlichen Daten wie Größe, Gewicht etc. ist das normalerweise der Fall, wenn das Modell richtig formuliert ist (d.h. wenn alle Faktoren, die einen Einfluss haben, enthalten sind und deren Beziehung - additiv, multitplikativ, exponentiell - korrekt dargestellt ist). Sind die Daten nicht kontinuierlich, können wir erstmal nicht von einer Normalverteilung der Residuen ausgehen. Häufig sind auch die Varianzen (= Streuung der Daten) nicht homogen über die Behandlungsgruppen. Das ist zum Beispiel bei Zähldaten der Fall, zumindest, wenn es sich um eher kleine Zahlen handelt (bis etwa 7).

4.2 Daten transformieren?

In so einem Fall wurde früher oft dazu geraten, die Daten zu transformieren, um sie ANOVA-fähig zu machen. Das birgt allerdings verschiedene Schwierigkeiten:

• die Interpretation der Ergebnisse wird schwieriger, weil die Ergebnisse und Signifikanzen nur für die transformierten Daten zutreffen. Damit verbunden ändert sich auch die Funktion des Zusammenhanges der erklärenden Parameter. Wenn die erklärenden Variablen vorher addiert wurden, um die abhängige Variable zu schätzen, müssen sie zum Beispiel nach einer Logarithmierung multipliziert werden.
• das Logarithmieren von Zähldaten, das häufig propagiert wird (wobei häufig eine 1 addiert wird, um Nullen zu vermeiden), führt häufig zu einer nicht gewollten Abnahme der Varianz der abhängigen Variable, die abhängig vom Mittelwert ist.
• in jedem Fall muss die Annahme der normalverteilten Residuen auch nach einer Transformation überprüft werden.

Die meisten Ratschläge, Daten zu transformieren, stammen aus der Zeit, als Varianzanalysen noch von Hand durchgeführt wurden und kompliziertere Modelle kaum zu rechnen waren. Zum Glück haben wir heute Computer und R, um auch viele nicht-normalverteilte Daten mittels einer ANOVA (Varianzanalyse) analysieren zu können. Der Trick hierbei ist, nicht die schon bekannten linearen Modelle (lm) zu verwenden, sondern generalisierte lineare Modelle (glm).

4.3 Generalisierte lineare Modelle

… können Poisson-verteilte Residuen (aus Zähldaten), binomial-verteilte Residuen (aus Anteilsdaten) und viele andere Residuenverteilungen beschreiben. Ich werde hier nicht auf die Charakteristika dieser Verteilungen eingehen, weil es für das Verständnis der Auswertung nicht essentiell ist. Wenn Sie Interesse haben, finden Sie aber gute Erklärungen über eine Internet-Suche.

Vergleich der Funktionsweisen eines linearen und eines generalisierten linearen Modells. Verändert nach einer Folie von Frank Schurr.

In einem normalen linearen Modell (oberer Teil der Abbildung) liefern die erklärenden Variablen (zum Beispiel ‘Düngung’)und die Modellgleichung direkt eine Schätzung der abhängigen Variable \(\mu\) (zum Beispiel ‘Ertrag). Die Abstände zu den tatsächlich beobachteten Daten \(y\) sind die Fehler oder Residuen und es wird angenommen, dass sie normalverteilt sind. Bei einem GLM wird zusätzlich eine ’link-function’ benötigt, die das Ergebnis der Modellgleichung auf den biologisch sinnvollen Bereich abbildet. So wird zum Beispiel verhindert, dass negative Zahlen geschätzt werden (ein Baum mit -5 Äpfeln). Zusätzlich ist die Annahme der Verteilung der Residuen/Fehler flexiber, wie oben bereits erwähnt.

Exkurs: Schätzung der Modell-Parameter
Grundsätzlich findet man die Modell-Parameter (bei einer Geradengleichung \(y = a * x + b\), wären das zum Beispiel \(a\) und \(b\)), für die die Daten am wahrscheinlichsten sind, indem man das Maximum-Likelihood Prinzip anwendet, sprich die deviance - also die Abweichung der geschätzten von den beobachteten Werten - minimiert. Bei normalverteilten Fehlern ist die deviance als die Summe der Fehlerquadrate RSS (residual sum of squares) definiert: \(\sum (\mu - y)^2\). Bei nicht-normalverteilten Fehlern funktioniert der Trick mit dem Minimieren der Summe der Fehlerquadrate nicht. Die deviance ist hier etwas komplizierter definiert. Für Poisson-verteilte Daten zum Beispiel \(2\sum(y_i log(y/\mu)-(y-\mu))\).

4.4 Beispiel: Tomaten

Leichter verständlich, wird es sicher, wenn wir uns ein konkretes Beispiel anschauen - hier können Sie einen Datensatz zu einem Tomatenversuch herunterladen: https://ecampus.uni-bonn.de/goto_ecampus_file_2786370_download.html Speichern Sie ihn in ihrem R-Projekt (wenn Sie eins angelegt haben) im Ordner ‘data’.

# Zuerst lesen wir den Tomaten-Datensatz ein
tomaten_daten <- read.csv("data/Tomaten.csv")

# Eine kurze Kontrolle, ob alles funktioniert hat
summary(tomaten_daten)

##     Datum             Zeitpunkt           Haus        Bedachung        
##  Length:865         Min.   : 1.000   Min.   :1.000   Length:865        
##  Class :character   1st Qu.: 3.000   1st Qu.:2.000   Class :character  
##  Mode  :character   Median : 5.000   Median :4.000   Mode  :character  
##                     Mean   : 5.089   Mean   :3.504                     
##                     3rd Qu.: 7.000   3rd Qu.:5.000                     
##                     Max.   :10.000   Max.   :6.000                     
##                                                                        
##   Pflanzen.Nr     Veredelung         Salzstress          Trieblänge   
##  Min.   : 1.00   Length:865         Length:865         Min.   :  0.0  
##  1st Qu.:25.00   Class :character   Class :character   1st Qu.: 80.0  
##  Median :49.00   Mode  :character   Mode  :character   Median :104.0  
##  Mean   :48.55                                         Mean   :106.4  
##  3rd Qu.:72.00                                         3rd Qu.:130.0  
##  Max.   :96.00                                         Max.   :196.0  
##                                                                       
##  Anzahl_Blätter  Anzahl_Blütenstände Anzahl_grüne_Früchte Anzahl_rote_Früchte
##  Min.   : 0.00   Min.   : 0.000      Min.   : 0.00        Min.   : 0.000     
##  1st Qu.:15.00   1st Qu.: 3.000      1st Qu.: 2.00        1st Qu.: 0.000     
##  Median :19.00   Median : 4.000      Median : 7.00        Median : 0.000     
##  Mean   :19.62   Mean   : 4.467      Mean   :11.14        Mean   : 1.029     
##  3rd Qu.:24.00   3rd Qu.: 6.000      3rd Qu.:17.00        3rd Qu.: 0.000     
##  Max.   :36.00   Max.   :36.000      Max.   :58.00        Max.   :15.000     
##                                                                              
##  Anzahl_Blütenendfäule Anzahl_geerntete_Früchte Gewicht_Trieb   Gewicht_Blätter
##  Min.   : 0.000        Min.   : 0.000           Min.   :  0.0   Min.   :   0   
##  1st Qu.: 0.000        1st Qu.: 0.000           1st Qu.:235.0   1st Qu.: 820   
##  Median : 1.000        Median : 0.000           Median :368.0   Median :1237   
##  Mean   : 3.042        Mean   : 0.341           Mean   :328.6   Mean   :1077   
##  3rd Qu.: 5.000        3rd Qu.: 0.000           3rd Qu.:428.0   3rd Qu.:1418   
##  Max.   :17.000        Max.   :10.000           Max.   :579.0   Max.   :1783   
##                                                 NA's   :768     NA's   :768    
##  Gewicht_grüne_Früchte Gewicht_rote_Früchte Gewicht_Blütenendfäule
##  Min.   :   0.0        Min.   :  0.0        Min.   :  0.00        
##  1st Qu.: 365.0        1st Qu.:  0.0        1st Qu.:  0.00        
##  Median : 688.0        Median :  0.0        Median :  0.00        
##  Mean   : 633.2        Mean   : 18.9        Mean   : 21.93        
##  3rd Qu.: 964.0        3rd Qu.:  0.0        3rd Qu.: 28.00        
##  Max.   :1222.0        Max.   :312.0        Max.   :221.00        
##  NA's   :768           NA's   :768          NA's   :768

# Es fällt auf, dass die beiden erklärenden Variablen 'Veredelung' und 'Salzstress', 
# die wir gleich verwenden wollen, as 'character' eingelesen worden sind. 
# Das kann bei den Abbildungen zu Problemen führen. Deshalb wandeln wir sie 
# zuerst in 'factors' um:
tomaten_daten$Veredelung <- as.factor(tomaten_daten$Veredelung)
tomaten_daten$Salzstress <-
  as.factor(tomaten_daten$Salzstress)

# Für die folgende Abbildung benötigen wir das package 'lattice' 
# (muss wie immer zuerst installiert und dann geladen werden).
library(lattice)

# Um die Daten ein bisschen kennenzulernen (und noch ein paar neue Arten von plots kennenzulernen), 
# bilden wir sie erstmal in einem xyplot ab. Insbesondere interessieren wir uns 
# jetzt für die Anzahl an Früchten, die Blütenendfäule haben in Abhängigkeit von 
# Salzstress und Veredelung.
# 
with(tomaten_daten, xyplot(Anzahl_Blütenendfäule ~ Salzstress, groups = Veredelung))

Hier sehen wir nicht besonders viel, außer dass die Varianz sehr hoch ist. Machen wir noch einen Interaktionsplot, in dem die Mittelwerte der Gruppen (mit/ohne Salzstress, mit/ohne Veredelung) berechnet und abgebildet werden.

# hier machen wir ein Interaktionsplot
with(tomaten_daten, interaction.plot(x.factor = Veredelung, 
   trace.factor = Salzstress, response = Anzahl_Blütenendfäule))

Jetzt erkennen wir, dass veredelte Tomaten mit Salzstress eine deutlich geringere Zahl von Früchten mit Blütenendfäule haben, als ohne Salzstress. Bei nicht-veredelten Pflanzen ist der Effekt genau umgekehrt, allerdings deutlich kleiner. Es sieht also so aus als gäbe es eine Interaktion zwischen den beiden erklärenden Variablen Salzstress und Veredelung: der Effekt der einen erklärenden Variablen auf die abhängige Variable hängt vom Wert der anderen erklärenden Variablen ab. Um diesen Effekte auf Signifikanz zu testen, fitten wir jetzt generalisierte lineare Modelle, die mit der Funktion glm() aufgerufen werden:

glm_model_1 <- glm(Anzahl_Blütenendfäule ~ Veredelung * Salzstress, 
    family = poisson, data = tomaten_daten)

Wie Sie sehen, funktioniert das ganz ähnlich wie mit der Funktion lm(). Es muss lediglich der zusätzliche Parameter ‘family’ angegeben werden. Hier wählen wir ‘poisson’, weil Zähldaten Poisson-verteilt sind.

4.5 Overdispersion

Aber Vorsicht! Bevor wir das Modell weiter analysieren können müssen wir uns noch das Verhältnis von Residual deviance (die Fehler, die es nach der Schätzung der Mittelwerte noch gibt) und den Freiheitsgraden/degrees of freedom anschauen: Ein GLM nimmt eine bestimmte Beziehung zwischen der Residual deviance und der Modellvorhersage an. Bei vielen agrarwissenschaftlichen Daten ist die Varianz aber größer als unter Poisson-Verteilung erwartet (z.B. weil nicht gemessene Variablen die abhängige Variable beeinflussen) -> das bezeichnet man als Overdispersion. Wenn die Residual deviance (zweitletzte Zeile in der summary) mehr als 1,5 mal höher ist als die Freiheitsgrade, liegt Overdispersion vor.

summary(glm_model_1)

## 
## Call:
## glm(formula = Anzahl_Blütenendfäule ~ Veredelung * Salzstress, 
##     family = poisson, data = tomaten_daten)
## 
## Coefficients:
##                               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                    1.01664    0.04093  24.840  < 2e-16 ***
## Veredelungohne                 0.08507    0.05675   1.499    0.134    
## Salzstressohne                 0.27847    0.05414   5.143 2.70e-07 ***
## Veredelungohne:Salzstressohne -0.37364    0.07854  -4.757 1.96e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 4293.7  on 864  degrees of freedom
## Residual deviance: 4256.1  on 861  degrees of freedom
## AIC: 5833.3
## 
## Number of Fisher Scoring iterations: 6

4256/861 ist deutlich größer als 1,5, wir haben in den Daten also tatsächlich Overdispersion.

Lösungen:

• wenn möglich, weitere erklärende Variablen mit in das Modell einbeziehen
• Wenn das nicht hilft/möglich ist: Verwendung von quasipoisson-Fehlern:

glm_model_1 <- glm(Anzahl_Blütenendfäule ~ Veredelung * Salzstress, 
      family = quasipoisson, data = tomaten_daten)

In diesem Modell wird nun zusätzlich geschätzt, wie hoch die Varianz in den einzelnen Behandlungsgruppen ist. Hiermit können wir weiterarbeiten: um die Interaktion zwischen Veredelung und Salzstress auf Signifikanz zu testen, formulieren wir ein einfacheres, genestetes Modell, in dem diese Interaktion fehlt (wenn sich ein Modell B von einem Modell A ableiten lässt, dass heißt, ein oder mehrere erklärende Variablen oder Interaktionen herausgenommen worden sind, sagt man, das Modell B ist in Modell A genestet). Dazu ersetzen wir das * (das für Interaktion steht) mit einem +. So haben beide erklärenden Variablen einen unabhängigen Einfluss auf die Anzahl der Früchte mit Blütenendfäule aber keine Interaktion mehr.

glm_model_2 <- glm(Anzahl_Blütenendfäule ~ Veredelung + Salzstress, 
        family = quasipoisson, data = tomaten_daten)

Mit der Funktion anova() vergleichen Sie die Modelle und testen, ob die Interkation signifikant ist. Zusätzlich müssen Sie hier die Statistik angeben, mit der die Modelle verglichen werden sollen, weil das bei glms mit eindeutig ist. Für quasipoisson-verteilte Residuen, wird der F-test benötigt.

anova(glm_model_1, glm_model_2, test = "F") 

## Analysis of Deviance Table
## 
## Model 1: Anzahl_Blütenendfäule ~ Veredelung * Salzstress
## Model 2: Anzahl_Blütenendfäule ~ Veredelung + Salzstress
##   Resid. Df Resid. Dev Df Deviance      F  Pr(>F)  
## 1       861     4256.1                             
## 2       862     4278.8 -1  -22.723 4.8598 0.02775 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Wie schon vermutet, ist die Interaktion signifikant. Wenn eine Interaktion signifikant ist, bleiben auch die beiden erklärenden Variablen, die interagieren auf jeden Fall im Modell. Die Varianzanalyse ist also jetzt abgeschlossen.

Wir könnten uns jetzt mit summary(*Modellname*) noch die geschätzten Parameter-Werte anschauen. Hier müssen Sie aber beachten, dass sie auf der link-Skala angegeben sind (siehe oben, die Parameter werden so angepasst, dass das Modell sinnvolle Ergebnisse liefert) und zurück-transformiert werden müssen, um den eigentlichen Wert zu erhalten. Am einfachsten ist es, die Funktion predict() zu verwenden: Als erstes Argument nennen wir den Namen des glm-Modells, dann erstellen wir einen Datensatz mit den Variablen und Werten, für die wir die Vorhersagen sehen möchten, also mit/ohne Salzstress und mit/ohne Veredelung in allen 4 möglichen Kombinationen. Als type wählen wir ‘response’ aus, weil wir die Schätzungen nicht auf der link-Skala sondern auf der response-Skala (das sind die zurücktransformierten, interpretierbaren Werte, siehe oben) haben wollen.

predict(glm_model_1, data.frame(Salzstress = c("mit", "mit", "ohne", "ohne"), 
                    Veredelung = c("mit", "ohne", "mit", "ohne")), type = "response")

##        1        2        3        4 
## 2.763889 3.009302 3.651376 2.736111

Der geschätzte Wert für die Anzahl an Früchten mit Blütenendfäule für Pflanzen mit Salzstress und mit Veredelung ist also 2,76, für Pflanzen mit Salzstress und ohne Veredelung 3.01 usw.

4.6 Anteils- und Boniturdaten auswerten

In dieser Woche geht es gleich um zwei Datentypen: Anteils- und Bonitur-Daten. Anteilsdaten können mit generalisierten linearen Modellen (GLMs) ausgewertet werden, die Sie beim letzten Mal schon kennengelernt haben - allerdings mit zwei kleinen Änderungen. Bonitur-Daten sind im Gartenbau ziemlich häufig, sind allerdings in Bezug auf die Auswertung ein eher undankbarer Datentyp: oft ist es schwierig bis unmöglich, sie durch Verteilungen wie der Normalverteilung zu charakterisieren. Deshalb wiedersetzen sie sich auch den Auswertungen, die auf solchen Verteilungen beruhen.

Am Ende dieses Kapitels wissen Sie

• wie die abhängige Variable in einem GLM zu Anteilsdaten generiert wird
• mit welchem Modell und welchem Test Sie Anteilsdaten auswerten, um signifikante Unterschiede zwischen den Behandlungen zu identifizieren
• wie man die etwas störrischen Bonitur-Daten zähmen, plotten und auswerten kann.

4.7 Anteilsdaten auswerten

Im Keimungsversuch mit Radis und Bohnen auf zwei verschiedenen Substraten haben Sie Anteilsdaten aufgenommen: wie viele der 80 ausgesäten Samen sind gekeimt? Die csv-Dateien finden Sie hier: https://ecampus.uni-bonn.de/goto_ecampus_file_2786369_download.html

Wie Sie sich schon denken können, erfüllt auch dieser Datentyp nicht die Annahmen, die für eine ‘normale’ ANOVa gegeben sein müssen. Zum Glück können wir wieder auf die generalisierten linearen Modelle zurückgreifen, allerdings mit zwei kleinen Änderungen. Aber fangen wir von vorne an.

Zuerst laden wir das Paket ‘tidyverse’, das uns zum Beispiel das Umstrukturieren von Datensätzen erleichtert, lesen die Daten ein (bei mir sind sie wieder im Ordner ‘data’ innerhalb des Ordners, in dem sich das r-Skript befindet gespeichert) und kontrollieren, ob alles geklappt hat:

library(tidyverse)

Ansaaten <- read.csv("data/Ansaaten.csv", header = T) 
head(Ansaaten)
str(Ansaaten)
summary(Ansaaten)

Um uns die Keimungsraten genauer anzuschauen, erstellen wir einen neuen Datensatz ‘Keimung’, indem wir mit dem pipe-Operator %>% aus dem Datensatz ‘Ansaaten’ die entsprechenden Spalten auswählen:

Keimung <- Ansaaten %>% select(Art, ID, Substrat, Ausfaelle, Angelaufen)
head(Keimung)

##     Art ID  Substrat Ausfaelle Angelaufen
## 1 Bohne  1 Presstopf         5         79
## 2 Bohne  2 Presstopf         3         81
## 3 Bohne  3 Presstopf         6         78
## 4 Bohne  4 Presstopf         6         78
## 5 Bohne  5 Presstopf        21         63
## 6 Bohne  6 Presstopf        23         61

Um einen Säulendiagramm zu erstellen, in dem die Anzahl der angelaufenen Samen und der Ausfälle gestapelt dargestellt werden, müssen wir den Datensatz umstrukturieren, um für jede Anzahl ‘Ausfaelle’ oder ‘Angelaufen’ eine eigene Spalte zu bekommen. Dazu nutzen wir die Funktion ‘reshape’. Wichtig sind hierbei die Parameter

• idvar = die Kombination aus Plot, ID und Behandlungslevels, die eine Zeile in der neuen Datenstruktur definiert
• timevar = Bezeichnung für neue Spalte
• times = das, was in diese Spalte eingetragen werden soll

Keimung_w <- reshape(Keimung, 
        direction = "long",
        varying = list(names(Keimung)[4:5]),
        v.names = "Anzahl",
        idvar = c("Art", "ID", "Substrat"),
        timevar = "Beobachtung",
        times = c("Ausfaelle", "Angelaufen"))
head(Keimung_w)

##                               Art ID  Substrat Beobachtung Anzahl
## Bohne.1.Presstopf.Ausfaelle Bohne  1 Presstopf   Ausfaelle      5
## Bohne.2.Presstopf.Ausfaelle Bohne  2 Presstopf   Ausfaelle      3
## Bohne.3.Presstopf.Ausfaelle Bohne  3 Presstopf   Ausfaelle      6
## Bohne.4.Presstopf.Ausfaelle Bohne  4 Presstopf   Ausfaelle      6
## Bohne.5.Presstopf.Ausfaelle Bohne  5 Presstopf   Ausfaelle     21
## Bohne.6.Presstopf.Ausfaelle Bohne  6 Presstopf   Ausfaelle     23

Jetzt können wir den Plot erstellen

ggplot(Keimung_w, aes(x=factor(ID), y=Anzahl, fill=Beobachtung)) +
  facet_grid(. ~ Art * Substrat) +
geom_bar(position="stack", stat="identity") +
  scale_fill_manual("legend", 
      values = c("Angelaufen" = "darkgreen", "Ausfaelle" = "lightgrey")) +
geom_col(position = position_stack(reverse = TRUE))

Wir sehen, dass die Erfolgsrate bei den Radis unabhängig vom Substrat sehr hoch ist. Bei den Bohnen gibt es mehr Ausfälle, insbesondere bei den auf Steinwolle ausgesäten.

Sind diese Unterschiede signifikant? Um das zu testen, können wir wieder die generalisierten linearen Modelle fitten und mittels der anova-Funktion vergleichen. Bei Anteilsdaten wie diesen erfolgt davor allerdings noch ein Vorbereitungsschritt: Wir generieren aus den beiden Spalten ‘Angelaufen’ und ‘Ausfaelle’ eine einzige abhängige Variable, die wir ‘Keimungsrate’ nennen. Es ist wichtig, die Rate nicht selber auszurechnen: ‘1 von 10 Samen gekeimt’ hat im Modell weniger Gewicht als ‘10 von 100’ Samen gekeimt. Stattdessen verwenden wir die Funktion ‘cbind()’ :

Ansaaten$Keimungsrate <- cbind(Ansaaten$Angelaufen, Ansaaten$Ausfaelle)

Dann können wir die Funktion glm() wie bekannt nutzen. Bei Anteilsdaten nutzen wir als family ‘binomial’

Modell1 <- glm(Keimungsrate ~ Substrat * Art, family = binomial, data = Ansaaten)
summary(Modell1)

## 
## Call:
## glm(formula = Keimungsrate ~ Substrat * Art, family = binomial, 
##     data = Ansaaten)
## 
## Coefficients:
##                              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                   1.61803    0.09285  17.427  < 2e-16 ***
## SubstratSteinwolle           -1.57041    0.11570 -13.574  < 2e-16 ***
## ArtRadies                     1.61079    0.20275   7.945 1.95e-15 ***
## SubstratSteinwolle:ArtRadies  1.33731    0.26856   4.980 6.37e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1266.76  on 39  degrees of freedom
## Residual deviance:  561.64  on 36  degrees of freedom
## AIC: 705.33
## 
## Number of Fisher Scoring iterations: 5

Die summary des Modells zeigt, dass die Daten overdispersed sind. Auch hier gibt es für diesen Fall die Möglichkeit als family ‘quasibinomial’ anzugeben:

Modell1 <- glm(Keimungsrate ~ Substrat * Art, family = quasibinomial, data = Ansaaten)

Zunächst lassen wir die Interaktion zwischen Substrat und Art weg, indem wir das Multiplikationszeichen durch ein Plus ersetzen

Modell2 <- glm(Keimungsrate ~ Substrat + Art, family = quasibinomial, data = Ansaaten)

Jetzt vergleichen wir die Modelle mittels anaova(). Als Test muss jetzt der F-Test angegeben werden

anova(Modell1, Modell2, test = "F")

## Analysis of Deviance Table
## 
## Model 1: Keimungsrate ~ Substrat * Art
## Model 2: Keimungsrate ~ Substrat + Art
##   Resid. Df Resid. Dev Df Deviance     F Pr(>F)
## 1        36     561.64                         
## 2        37     584.97 -1  -23.326 1.603 0.2136

Wie wir sehen, ist die Interaktion nicht signifikant (p > 0.05) Also vereinfachen wir das Modell weiter und nehmen jetzt das Substrat als erklärende Variable raus

Modell3 <- glm(Keimungsrate ~ Art, family = quasibinomial, data = Ansaaten)
# Wieder vergleichen:
anova(Modell2, Modell3, test = "F")

## Analysis of Deviance Table
## 
## Model 1: Keimungsrate ~ Substrat + Art
## Model 2: Keimungsrate ~ Art
##   Resid. Df Resid. Dev Df Deviance      F   Pr(>F)   
## 1        37     584.97                               
## 2        38     767.96 -1  -182.99 11.644 0.001573 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Das Substrat hat einen hochsignifikanten Effekt auf die Keimungsrate. Als letztes testen wir noch die Art

Modell4 <- glm(Keimungsrate ~ Substrat, family = quasibinomial, data = Ansaaten)

Jetzt müssen wir Modell4 natürlich gegen Modell2 testen. Die Modelle müssen immer ‘genestet’ sein, sonst kann kein Signifikanzniveau berechnet werden. Genestet bedeutet, dass ein Modell eine einfachere Variante (eine oder mehrere Variable/n oder Interaktion/en weniger) des anderen Modells ist.

anova(Modell2, Modell4, test = "F")

## Analysis of Deviance Table
## 
## Model 1: Keimungsrate ~ Substrat + Art
## Model 2: Keimungsrate ~ Substrat
##   Resid. Df Resid. Dev Df Deviance      F    Pr(>F)    
## 1        37     584.97                                 
## 2        38    1108.34 -1  -523.38 33.303 1.283e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Auch die Art hat einen hochsignifikanten Einfluss auf die Keimungsrate. Damit ist die Analyse der Keimungsraten abgeschlossen.

4.8 Bonitur-Daten

Bonitur-Daten sind recht speziell, weil es keine metrischen Daten sind, sondern Ordinal-Daten und darüber hinaus auch einigermaßen subjektiv. Trotzdem sind hierfür Methoden zur Auswertung entwickelt worden. Sie stammen aus der Psychologie, in der häufig Untersuchungen stattfinden, bei denen Personen zum Beispiel ihre Zufriedenheit auf einer Skala von 1 bis 5 bewerten sollen.

Wir schauen uns dazu exemplarisch die Bonitur-Daten zur Entwicklung der Bohnen an. Dazu filtern wir den Datensatz Ansaaten, indem wir nur die Zeilen nehmen, in denen das Argument Art == “Bohne” zutrifft (das doppelte Gleichzeichen wird immer benutzt, wenn eine Abfrage gemacht wird, ein einfaches Gleichzeichen ist in R eine Zuweisung also synonym mit dem Pfeil <- ). Mit dem Pipe-Operator %>% verbinden wir diese Auswahl mit der Auswahl der Spalten, die wir mit der Funktion select() machen. Hier können wieder einfach die header (=Spaltenüberschriften) der benötigten Spalten angegeben werden. Vorher laden wir noch zwei Pakete, die für die Auswertung von Ordinal-Daten entwickelt wurden

library(likert)

## Loading required package: xtable

## 
## Attaching package: 'likert'

## The following object is masked from 'package:car':
## 
##     recode

## The following object is masked from 'package:dplyr':
## 
##     recode

library(ordinal)

## 
## Attaching package: 'ordinal'

## The following object is masked from 'package:dplyr':
## 
##     slice

Bohne_Entw_l <- filter(Ansaaten, Art == "Bohne") %>% select(Substrat, ID, Entwicklung)
head(Bohne_Entw_l)

##    Substrat ID Entwicklung
## 1 Presstopf  1           8
## 2 Presstopf  2           5
## 3 Presstopf  3           8
## 4 Presstopf  4           8
## 5 Presstopf  5           4
## 6 Presstopf  6           5

Als erstes sollten die Daten wieder geplottet werden. Bei Bonitur-Daten eignen sich Histogramme am besten. Um das zu tun, müssen die Bonitur-Noten in Faktoren umgewandelt werden. Für die weitere Analyse geben wir auch direkt an, wieviele mögliche Levels (verschiedene Noten) es gibt und dass sie eine aussagekräftige Reihenfolge haben (ordered = TRUE):

Bohne_Entw_l$Entwicklung <- factor(Bohne_Entw_l$Entwicklung, ordered = TRUE, levels = 1:9)

Jetzt können die Histogramme erstellt werden:

ggplot(Bohne_Entw_l, aes(Entwicklung)) +
  geom_bar() +
  facet_wrap(~Substrat, ncol=1) +
  xlab("Bonitur Note") +
  ylab("Anzahl")

Das package ‘likert’ stellt einige Plot-Funktionen zur Verfügung, die speziell für Bonitur- und andere Typen von likert-Daten geeignet sind. Dazu brauchen wir die Daten allerdings im ‘wide’ Format. Im Moment haben wir die Daten im ‘long’ Format: es gibt eine eigene Spalte für das Substrat, sodass ‘Presstopf’ und ‘Steinwolle’ je 10 Mal genannt werden. Wir nutzen wieder die Funktion ‘reshape’ wie oben, jetzt allerdings in umgekehrter Richtung, vom ‘long’ zum ‘wide’ Format:

Bohne_Entw_w <- reshape(Bohne_Entw_l, v.names="Entwicklung", idvar = "ID", 
    timevar = "Substrat", direction="wide")
head(Bohne_Entw_w)

##   ID Entwicklung.Presstopf Entwicklung.Steinwolle
## 1  1                     8                      5
## 2  2                     5                      5
## 3  3                     8                      5
## 4  4                     8                      3
## 5  5                     4                      2
## 6  6                     5                      1

‘idvar’ sind jetzt also die Zeilenbezeichnungen, ‘timevar’ die Spaltenbezeichnungen, ‘v.names’ die eigentlichen Werte.

Als nächstes müssen wir ein sogenanntes ‘likert-Objekt’ aus den beiden Spalten mit den Bonitur-Noten erstellen, dann kann die plot-Funktion des likert-packages genutzt werden:

Bohne_Entw_likert = likert(Bohne_Entw_w[2:3])

plot(Bohne_Entw_likert)

Dies ist eine andere - etwas kompaktere - Darstellung der gleichen Daten. Die %-Werte geben an, wieviel % im schlechten Bereich (1-3), im mittleren Bereich (4-6) und im guten Bereich (7 - 9) waren.

Bonitur-Daten sind Ordinal-Daten und haben damit andere Eigenschaften als metrische Daten. Wie oben erwähnt, erfüllen sie häufig nicht die Annahmen, die für eine normale Varianzanalyse zutreffen müssen. Wir greifen hier auf ‘cumulative link models’ (CLM) aus der library ‘ordered’ zurück. CLMs sind im Prinzip das gleiche wie proportional odds model, falls Ihnen dieser Begriff mal über den Weg läuft. Ähnlich wie in den vorherigen beiden Kapiteln formulieren wir wieder ein komplexeres Modell mit Substrat als erklärender Variable und das Null-Modell, indem angenommen wird, dass alle Bohnen sich ungeachtet der Behandlung gleich Entwickeln:

Model1 <- clm(Entwicklung ~ Substrat, data = Bohne_Entw_l, threshold = "equidistant")

Model2 <- clm(Entwicklung ~ 1, data = Bohne_Entw_l, threshold = "equidistant")

Als threshhold geben wir ‘equidistant’ an, weil wir davon ausgehen, dass zum Beispiel der Unterschied in der Entwicklung zwischen den Bonitur-Noten 2 und 3 genauso groß ist, wie der Unterschied zwischen den Noten 3 und 4 usw.

Auch für CLMs gibt es die Funktion anova(). Hier wird allerdings kein t- oder F-Test durchgeführt, sondern ein Chi-square Test, der uns einen p-Wert liefert:

anova(Model1, Model2)

## Likelihood ratio tests of cumulative link models:
##  
##        formula:               link: threshold: 
## Model2 Entwicklung ~ 1        logit equidistant
## Model1 Entwicklung ~ Substrat logit equidistant
## 
##        no.par    AIC  logLik LR.stat df Pr(>Chisq)   
## Model2      2 78.001 -37.000                         
## Model1      3 73.330 -33.665  6.6711  1   0.009799 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Wie sie sehen, ist der p-Wert etwa 0.01, das heißt, das Substrat hat einen signifikanten Effekt auf die Entwicklung der Bohnen. Hiermit ist auch eine grundlegende Analyse der Bonitur-Daten abgeschlossen.

Ein interessanter Link zum Thema, inklusive Auswertungsmöglichkeiten mit R: ‘Do not use averages with Likert scale data’

5.1 Zeitreihen plotten

Zunächst laden wir die library ‘ggplot2’, lesen die Wachstumsdaten aus der Datei ‘Braeburn.csv’ ein (im Order ‘data’ auf eCampus) und kontrollieren den Import. ID und Zeitpunkt wandeln wir in Faktoren um. Wie Sie sehen, haben wir Daten von 4 Früchten, 2 von Ästen mit niedrigen Behangdichten und 2 von Ästen mit hohen Behangdichten, deren Umfang stündlich gemessen wurde.

library(ggplot2)

Braeburn <- read.csv("data/Braeburn.csv", header = TRUE)
summary(Braeburn)

##     Datum               Umfang            ID          Behang         
##  Length:2508        Min.   : 1434   Min.   :1.00   Length:2508       
##  Class :character   1st Qu.: 5049   1st Qu.:1.75   Class :character  
##  Mode  :character   Median : 7294   Median :2.50   Mode  :character  
##                     Mean   : 7443   Mean   :2.50                     
##                     3rd Qu.: 9583   3rd Qu.:3.25                     
##                     Max.   :14609   Max.   :4.00                     
##    Zeitpunkt  
##  Min.   :  1  
##  1st Qu.:157  
##  Median :314  
##  Mean   :314  
##  3rd Qu.:471  
##  Max.   :627

Braeburn$ID <- as.factor(Braeburn$ID)
Braeburn$Zeitpunkt <- as.factor(Braeburn$Zeitpunkt)

Um das Wachstum dieser 4 Früchte zu visualisieren, nutzen wir die ggplot Funktion und geben bei den Aestetics (ein Parameter der Funktion) den Zeitpunkt (x-Achse), den Umfang dividiert durch 1000 um von Mikro- auf Millimeter zu kommen (y-Achse) und den Behang (niedrig/hoch) für die Gruppierung durch Farbe an:

ggplot(Braeburn, aes(Zeitpunkt, Umfang/1000, color=Behang)) + 
   geom_point() +   labs(y="Umfang (mm)", x="Zeit") +
   theme_bw()

Im besten Fall, insbesondere wenn wir Unterschiede im Wachstum durch verschiedene Behandlung statistisch nachweisen wollen, haben wir nicht nur 2 Wiederholungen wie hier sondern deutlich mehr. Dies ist im Datensatz ‘Hühner’ der Fall, in dem das Gewicht von Junghennen, die mit unterschiedlichem Futter aufgezogen wurden, über 21 Tage protokolliert wurde.

Auch diesen Datensatz lesen wir ein (auch Ordner ‘data) und wandeln ’Futter’ und ‘ID’ in Faktoren um:

Hühner <- read.csv("data/Huehner.csv", header = TRUE)
head(Hühner)

##   X Gewicht Tag ID Futter
## 1 1      42   0  1      1
## 2 2      51   2  1      1
## 3 3      59   4  1      1
## 4 4      64   6  1      1
## 5 5      76   8  1      1
## 6 6      93  10  1      1

summary(Hühner)

##        X            Gewicht           Tag              ID       
##  Min.   :  1.0   Min.   : 35.0   Min.   : 0.00   Min.   : 1.00  
##  1st Qu.:145.2   1st Qu.: 63.0   1st Qu.: 4.00   1st Qu.:13.00  
##  Median :289.5   Median :103.0   Median :10.00   Median :26.00  
##  Mean   :289.5   Mean   :121.8   Mean   :10.72   Mean   :25.75  
##  3rd Qu.:433.8   3rd Qu.:163.8   3rd Qu.:16.00   3rd Qu.:38.00  
##  Max.   :578.0   Max.   :373.0   Max.   :21.00   Max.   :50.00  
##      Futter     
##  Min.   :1.000  
##  1st Qu.:1.000  
##  Median :2.000  
##  Mean   :2.235  
##  3rd Qu.:3.000  
##  Max.   :4.000

Hühner$Futter <- as.factor(Hühner$Futter)
Hühner$ID <- as.factor(Hühner$ID)

Das entsprechende Diagramm sieht folgendermaßen aus:

ggplot(Hühner, aes(Tag, Gewicht, color=Futter)) + 
   geom_point() +   labs(y="Gewicht (g)", x="Tage seit Geburt") +
   theme_bw()

Für eine bessere Übersicht wäre es wünschenswert, nicht alle Datenpunkte zu plotten¸ sondern den Mittelwert des Gewichts der Hühner einer Behandlung und ein Maß entweder für die Varianz in den Daten oder die Genauigkeit des Mittelwertes. Mit der ggplot-Funktion ‘stat_summary’ können wir genau das machen: der Wert ‘mean_se’ im Parameter ‘fun.data’ bewirkt, dass Mittelwert und Standardfehler des Mittelwertes berechnet und geplottet werden:

ggplot(Hühner, aes(Tag, Gewicht, color=Futter)) + 
  stat_summary(fun.data=mean_se, geom="pointrange") + 
  labs(y="Gewicht (g)", x="Zeit (Tage seit Geburt)")

Wir sehen in diesem Plot einen deutlichen Trend, nämlich dass das Futter des Typs 3 zu schnellerer Gewichtszunahme führt, als die übrigen. Die nächsten beiden Abschnitte erklären, wie wir testen, ob dieser Trend tatsächlich signifikant ist.

5.2 Zufällige Effekte

Im Prinzip sollte man meinen, dass sich Gewichts- oder Wachstumsdaten gut eignen sollten, um eine Varianzanalyse durchführen zu können: sie sind kontinuierlich und wir erwarten im Prinzip eine Normalverteilung der Residuen. Allerdings ist eine grundlegende Annahme verletzt: die Unabhängigkeit der Daten. Wenn wir immer wieder den gleichen Apfel oder das gleiche Huhn messen, können wir davon ausgehen, dass die Werte einander ähnlicher sind, als wenn wir immer andere Indiviuen nehmen würden. Wir haben also keine echten unabhängigen Wiederholungen sondern Pseudoreplikationen. Andere Beispiele die zu Pseudoreplikationen führen sind Versuchsblöcke, Standorte, unterschiedliche Probennehmer, etc.

Wie können solche Faktoren adäquat im Modell berücksichtigt werden? Der klassische Ansatz war lange, sie als ‘normale’ Variablen, also als festen Effekt mit ins Modell aufzunehmen. Das verbraucht aber viele Freiheitsgrade, weil für jeden Level des Zufallseffekts (also jedes Huhn) ein Parameter geschätzt wird (z.B. das mittlere Gewicht jedes Huhns, siehe Abbildung unten). Hier gibt es eine recht gute Erklärung für Freiheitsgrade: https://welt-der-bwl.de/Freiheitsgrade. Je weniger Freiheitsgrade es gibt, desto schwieriger ist es in einer Varianzanalyse, tatsächlich existierende Unterschiede in den Behandlungsgruppen auch als solche zu erkennen (Fähigkeit, Unterschiede zu erkennen = Power der Analyse). Deshalb ist der bessere Ansatz, die Variable als Zufallseffekt einzubeziehen. Das beruht auf der Annahme, dass die Basis-Parameter für jedes Level dieser Variable (jedes Huhn) aus einer gemeinsamen Verteilung stammen. Geschätzt werden muss für das Modell nur die Breite, also Varianz dieser Verteilung. Auf diese Weise wird nur 1 Freiheitsgrad verwendet und die Power der Analyse bleibt höher.

Linke Seite: Klassischer Ansatz - für jedes Huhn wird ein Parameter für den Achsenabschnitt Gewicht geschätzt. Rechte Seite: Nur die breite der Verteilung, aus der die Parameter stammen, wird geschätzt, es wird also nur 1 Freiheitsgrad ‘verbraucht’. Verändert nach einer Folie von Frank Schurr

Feste Effekte

• die geschätzten Paremeter sind von Interesse, es sind zum Beispiel die unterschiedlichen Behandlungslevels (Behangdichte, Futter, …)
• die Namen der Levels haben eine (agrarwissenschaftliche) Bedeutung

Zufällige Effekte

• die geschätzten Parameter interessieren uns nicht
• die Namen der Level haben keine agrarwissenschafliche Bedeutung
• nur einbezogen, um Pseudoreplikation zu vermeiden (ID des Huhn, des Apfels, Block, Ort, Probennehmer)

5.3 Zeitreihen analysieren

Um zufällige Faktoren zu berücksichtigen, benötigen wir gemischte Modelle.

Vergleich der Struktur eines linearen Modells mit einem gemischten linearen Modell, verändert nach einer Folie von Frank Schurr

Ein Paket, das Funktionen für gemischte Modelle zur Verfügung stellt ist ‘lme4’. Dieses Paket laden wir jetzt (nachdem Sie es installiert haben).

library(lme4)

## Loading required package: Matrix

## 
## Attaching package: 'Matrix'

## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack

Die Funktion lmer() fitted lineare gemischte Modelle. Die abhängige Variable und die festen erklärenden Variablen werden wie bekannt in Zusammenhang gebracht: In unserem Beispiel steht das Gewicht links von der Tilde und Tag und Futter werden mit einem Multiplikationszeichen verbunden, um auszudrücken, dass in dem Modell beide Faktoren und ihre Interaktion einen Einfluss auf das Gewicht haben können. Dahinter kommen in runden Klammern die Zufallseffekte: Hier gehen wir davon aus, dass es einen Zufallseffekt geben könnte (individuelle Unterschiede im Gewicht der Hühner), der sich auf den y-Achsenabschnitt der Zeitreihen auswirkt. Das wird ausgedrückt, indem die ID hinter den senkrechten Strich geschrieben wird. Um die Modelle im Anschluss vergleichen zu können, müssen wir noch die Standardmethode ‘REML’ auf ‘FALSE’ setzen.

model1 <- lmer(Gewicht ~ Tag * Futter + (1| ID), data=Hühner, REML = FALSE)

Danach geht es wie gewohnt weiter: wir vereinfachen das Modell, zuerst durch Ersetzen des Multiplikationszeichens mit einem Additionszeichen, womit wir die Interaktion zwischen Tag und Futter rauswerfen. Ist die Interaktion nicht signifikant, nehmen wir als nächstes das Futter ganz raus. Der Zeitpunkt (=Tag) sollte bei solchen Zeitreihen-Analysen allerdings immer im Modell bleiben, es sei denn, wir gehen davon aus, dass es möglicherweise gar keine Veränderung der abhängigen Variable über die Zeit gibt.

model2 <- lmer(Gewicht ~ Tag + Futter + (1 | ID), data=Hühner, REML = F)
model3 <- lmer(Gewicht ~ Tag + (1 | ID), data=Hühner, REML = F)

Dann kann man die Modelle mittels anova() vergleichen. Wie Sie sehen, können auch mehrere Modelle gleichzeitig verglichen werden.

anova(model1, model2, model3)

## Data: Hühner
## Models:
## model3: Gewicht ~ Tag + (1 | ID)
## model2: Gewicht ~ Tag + Futter + (1 | ID)
## model1: Gewicht ~ Tag * Futter + (1 | ID)
##        npar    AIC    BIC  logLik -2*log(L)   Chisq Df Pr(>Chisq)    
## model3    4 5630.3 5647.8 -2811.2    5622.3                          
## model2    7 5619.2 5649.7 -2802.6    5605.2  17.143  3  0.0006603 ***
## model1   10 5508.0 5551.6 -2744.0    5488.0 117.184  3  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Sie sehen, dass sowohl der Einfluss des Futters auf den Gewichtsverlauf als auch die Interaktion zwischen Futter und Tag hochsignifikant sind. Die Vorhersagen des besten Modells (also model1) können wir nun noch zu dem Plot hinzufügen:

last_plot() + stat_summary(aes(y=fitted(model1)), fun=mean, geom="line")

5.4 Post-hoc Tests

Nachdem wir in den letzten Kapiteln Varianzanalysen mit verschiedenen Datentypen besprochen haben, können Sie testen, ob eine erklärende Variable “signifikant” ist, also ihr Einbezug in ein Modell die übrigbleibenden Fehler (=Residuen) signifikant reduziert. Was Sie noch nicht wissen ist, welche der Level (=Behandlungsgruppen) der erklärenden Variable signifikant unterschiedlich sind. Zum Beispiel haben Sie beim Ansaaten-Versuch den Effekt von 4 verschiedene Substraten auf die Keimungsrate von 2 Arten getestet. Mit einer ANOVA finden Sie heraus, ob das Substrat grundsätzlich einen signifikaten Effekt auf die Keimungsrate hat. Bei WELCHEM Substrat oder welchen Substraten aber die Keimungsrate signifikant höher ist oder sind als bei den anderen, wissen Sie noch nicht. Um das herauszufinden, benutzen wir post-hoc Tests, also Tests, die nach der ANOVA durchgeführt werden.

Am Ende dieses Kapitels wissen Sie

• um die Problematik der multiplen post-hoc Tests
• wie Sie solche Tests sinnvoll durchführen können

5.5 Post-hoc Tests - die Problematik

Natürlich können wir jetzt jedes Substrat nochmal einzeln mit einer ANOVA gegen jedes andere Substrat testen und schauen, ob das Ergebnis siginifkant ist. Allerdings gibt es dabei ein Problem: jeder Test hat per Definition eine Irrtumswahrscheinlichkeit von 5% (wir nennen eine erklärende Variable signifikant, wenn der p-Wert, also die Wahrscheinlichkeit solche oder noch extremere Daten zu finden, obwohl die Null-Hypothese wahr ist, kleiner als 0.05 = 5% ist). Je mehr Vergleiche wir machen, desto höher ist also auch die Wahrscheinlichkeit, dass wir in der gesamten Auswertung eine Signifikanz finden, obwohl gar kein wahrer Effekt vorhanden ist - das Problem der multiplen Tests.

source: https://xkcd.com/882/

Einige Statistiker sagen deshalb, dass post-hoc Tests grundsätzlich vermieden werden sollten und auch nicht nötig sind, wenn das experimentelle Design des zugrundeliegenden Versuchs nur gut genug angelegt wurde. Die agrarwissenschaftlichen Realität sieht aber oft anders aus: es wird häufig eine Vielzahl von Behandlungslevels gegeneinander getestet und es interessiert uns nicht nur, OB die Wahl des Substrats einen Effekt auf die Keimungsrate einer Art hat, sondern auch bei welchem davon die Keimungsrate signifikant höher ist. Wahr ist aber, dass man gut durchdenken sollte, wann post-hoc Tests wirklich gebraucht werden und wie genau man sie durchführt.

5.6 Die Herangehensweise

Die generelle Herangehensweise ist, dass man für die Anzahl an Vergleichen, die man macht, korrigiert, indem die einzelnen Signifikanz-Niveaus so angepasst werden, dass man INSGESAMT auf eine Irrtumswahrscheinlichkeit von 5% kommt. Gängige Ansätze sind die Bonferroni-Korrektur und der Tukey’s HSD Test. Hier nutzen wir eine angepasste Bonferroni Korrektur, die sich ‘Holm’ nennt. Wer sich für die methodischen Details interessiert findet hier eine recht gute Erklärung.

Bitte laden Sie die Daten zum Ansaaten-Versuch von diesem Jahr herunter: https://ecampus.uni-bonn.de/goto_ecampus_file_2859385_download.html
und installieren Sie das R-Paket ‘multcomp’.

# Laden der benötigten Pakete
library(multcomp)

## Loading required package: mvtnorm

## Loading required package: survival

## Loading required package: TH.data

## Loading required package: MASS

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:dplyr':
## 
##     select

## 
## Attaching package: 'TH.data'

## The following object is masked from 'package:MASS':
## 
##     geyser

library(tidyverse)

# Einlesen der Daten und Umwandeln der erklärenden Variablen in Faktoren
ansaaten <- read.csv("data/Ansaaten_neu.csv", header = TRUE)
ansaaten$Art <- as.factor(ansaaten$Art)
ansaaten$Substrat <- as.factor(ansaaten$Substrat)

# Kombinieren der beiden Spalten mit Anteilsdaten
ansaaten$keimung <- cbind(ansaaten$Angelaufen, ansaaten$Ausfälle)

# Umformen des Datensatzes, um einen Plot zu machen
Keimung_w <- reshape(ansaaten, 
                     direction = "long",
                     varying = list(names(ansaaten)[4:5]),
                     v.names = "Anzahl",
                     idvar = c("Art", "Substrat", "Wdh"),
                     timevar = "Beobachtung",
                     times = c("Ausfaelle", "Angelaufen"))
head(Keimung_w)

##                                Art   Substrat Wdh keimung.1 keimung.2
## Bohne.Jiffy.1.Ausfaelle      Bohne      Jiffy   1        63        21
## Bohne.Jiffy.2.Ausfaelle      Bohne      Jiffy   2        34        50
## Bohne.Presstopf.1.Ausfaelle  Bohne  Presstopf   1        76         8
## Bohne.Presstopf.2.Ausfaelle  Bohne  Presstopf   2        76         8
## Bohne.Steinwolle.1.Ausfaelle Bohne Steinwolle   1        64        20
## Bohne.Steinwolle.2.Ausfaelle Bohne Steinwolle   2        70        14
##                              Beobachtung Anzahl
## Bohne.Jiffy.1.Ausfaelle        Ausfaelle     63
## Bohne.Jiffy.2.Ausfaelle        Ausfaelle     34
## Bohne.Presstopf.1.Ausfaelle    Ausfaelle     76
## Bohne.Presstopf.2.Ausfaelle    Ausfaelle     76
## Bohne.Steinwolle.1.Ausfaelle   Ausfaelle     64
## Bohne.Steinwolle.2.Ausfaelle   Ausfaelle     70

# Plot
ggplot(Keimung_w, aes(x=factor(Wdh), y=Anzahl, fill=Beobachtung)) +
  facet_grid(. ~ Art * Substrat) +
  geom_bar(position="stack", stat="identity") +
  scale_fill_manual("legend", 
                    values = c("Angelaufen" = "darkgreen", "Ausfaelle" = "lightgrey")) +
  geom_col(position = position_stack(reverse = TRUE))

# Generalisiertes lineares Modell unserer Hypothese: 
# Substrat und Art haben einen Effekt auf die Keimungsrate
model1 <- glm(keimung ~ Substrat * Art, family=quasibinomial, data = ansaaten)

# Genestetes Modell ohne den Effekt von 'Substrat', 
# um die Signifikanz dieser Variable testen zu können
model2 <- glm(keimung ~ Art, family = quasibinomial, data = ansaaten)

# Varianzanalyse der beiden Modelle
anova(model1, model2, test = "F")

## Analysis of Deviance Table
## 
## Model 1: keimung ~ Substrat * Art
## Model 2: keimung ~ Art
##   Resid. Df Resid. Dev Df Deviance      F    Pr(>F)    
## 1         8      31.78                                 
## 2        14     326.94 -6  -295.16 12.977 0.0009685 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Mit der Funktion interaction() generieren wir jetzt eine neue Spalte ‘gruppe’, in der die erklärenden Variablen kombiniert werden. Mit diesen Behandlungsgruppen fitten wir ein neues glm, um sie dann in einem nächsten Schritt in einer post-hoc Analyse gegeneinander testen zu können.

ansaaten$gruppe<-interaction(ansaaten$Art, ansaaten$Substrat)
head(ansaaten)

##     Art   Substrat Wdh Angelaufen Ausfälle keimung.1 keimung.2           gruppe
## 1 Bohne      Jiffy   1         63       21        63        21      Bohne.Jiffy
## 2 Bohne      Jiffy   2         34       50        34        50      Bohne.Jiffy
## 3 Bohne  Presstopf   1         76        8        76         8  Bohne.Presstopf
## 4 Bohne  Presstopf   2         76        8        76         8  Bohne.Presstopf
## 5 Bohne Steinwolle   1         64       20        64        20 Bohne.Steinwolle
## 6 Bohne Steinwolle   2         70       14        70        14 Bohne.Steinwolle

model_posthoc <- glm(keimung ~ gruppe, family=quasibinomial, data = ansaaten)

Bevor wir das tun, sollten wir überlegen, welche Gruppen wir tatsächlich vergleichen wollen. Die Abbildung oben legt nahe, dass Jiffy stripes zu einer höheren Keimungsrate führen als Steinwolle und Sand und wir wollen vielleicht testen, ob diese Unterschiede signifikant sind. In der Funktion glht , die den post-hoc Test für uns durchführt können wir festlegen, welche Gruppen wir gegeneinander testen wollen:

summary(glht(model_posthoc, 
             linfct = mcp(gruppe =
#Ist die Differenz zwischen diesen Gruppen ungleich 0?
            c("(Bohne.Jiffy)-(Bohne.Sand)=0",
              "(Bohne.Steinwolle)-(Bohne.Sand)=0"))),test = adjusted("holm"))

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: User-defined Contrasts
## 
## 
## Fit: glm(formula = keimung ~ gruppe, family = quasibinomial, data = ansaaten)
## 
## Linear Hypotheses:
##                                        Estimate Std. Error z value Pr(>|z|)  
## (Bohne.Jiffy) - (Bohne.Sand) == 0       -1.4319     0.5202  -2.753   0.0118 *
## (Bohne.Steinwolle) - (Bohne.Sand) == 0  -0.3725     0.5639  -0.661   0.5089  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- holm method)

Jetzt sehen Sie, dass der Unterschied zwischen Jiffy stripes und Sand bei der Bohne signifikant ist, der Unterschied Steinwolle - Sand aber nicht.

Signifikante Unterschiede, die zwischen den Behandlungsgruppen gefunden wurden, werden in Abbildungen häufig durch verschiedene Buchstaben gekennzeichnet. Zum Beispiel so (hier liegen andere Daten zugrunde):

https://r-graph-gallery.com/84-tukey-test.html

Allerdings bietet R keine guten Standardlösungen dafür an. Wenn Sie den TukeyHSD Test durchführen, bietet das Packet ‘multcompView’ eine Möglichkeit. R-Code um die Abbildung oben zu erstellen finden Sie hier.

6.1 Der Satz von Bayes

Keine Einführung in die Bayesische Statistik kommt ohne den Satz von Bayes aus: Thomas Bayes war ein englischer Geistlicher, Philosoph und Statistiker von dem 1763 - posthum - ein Manuskript mit dem Titel ‘An essay towards solving a problem in the doctrine of chances’ veröffentlicht wurde. In diesem beschreibt er, wie basierend auf vorhandenem Wissen und zusätzlichen Beobachtungen X die Wahrscheinlichkeit dafür berechnet werden kann, dass eine Hypothese H zutrifft:

\[P(H|X) = \frac {P(X|H)×P(H)}{P(X)}\]

In Worten ausgedrückt: Die Wahrscheinlichkeit (P = probability) dass die Hypothese zutrifft, nachdem wir die Daten betrachtet haben (P(H|X), a-posteriori Wissen), ist gleich der Wahrscheinlichkeit, die Daten so zu beobachten, angenommen die Hypothese trifft zu (P(X|H), Likelihood), mal der Wahrscheinlichkeit, dass die Hypothese zutrifft, bevor wir die Daten betrachtet haben (P(H), a-priori Wissen oder Prior), geteilt durch die Wahrscheinlichkeit die Daten so zu beobachten ohne irgendeine Annahme zur Hypothese (P(X), Normalisierungskonstante).
Wichtig ist noch anzumerken, dass Hypothesen dabei als Parameterwerte ausgedrückt werden. Ist unsere Hypothese also, dass unter Behandlung A der Mittelwert der abhängigen Variable größer ist als unter Behandlung B, berechnen wir die Wahrscheinlichkeit für den Parameter Mittelwert(A) - Mittelwert(B) und bekommen als Ergebnis, mit welcher Wahrscheinlichkeit er > 0 ist.

6.2 Das Prinzip des Bayes-Statistik

In der Bayes-Statistik wird diese Idee verwendet, um bereits vorhandenes Wissen - das A-priori-Wissen (P(HP)) - mit den Information, die in den neu erhobenen Daten X stecken zu kombinieren und so das ‘geupdatete- A-posteriori-Wissen zu generieren. Wir sehen außerdem, dass wir als Ergebnis die Wahrscheinlichkeit unserer Hypothese bekommen. Das ist wesentlich leichter bekömmlich und näher an unserem ’normalen’ Denken als die Interpretation eines p-Wertes: ‘Die Wahrscheinlichkeit solche oder extremere Daten zu finden, wenn die Null-Hypothese wahr ist’.

Exkurs: Wahrscheinlichkeitsverteilungen
Um Wahrscheinlichkeiten darzustellen, zum Beispiel das a-priori Wissen, werden Wahrscheinlichkeitsverteilungen verwendet. In der foldenden Abbildung sehen Sie eine solche Wahrscheinlichkeitsverteilung: Je höher die Dichte, desto wahrscheinlicher der Parameterwert, der auf der x-Achse abgetragen ist. Die Gesamtfläche unter der Kurve beträgt immer 1 (genau irgendein Parameterwert trifft immer zu).

Wie berechnen wir also diese Wahrscheinlichkeit? Leider ist das meistens nicht ganz so trivial, wie es aussieht. Zwar können wir meistens die Likelihood berechnen (das wird auch in der frequentistischen Statistik zur Bestimmung der Parameterwerte unserer Modelle gemacht) und unser a-priori Wissen durch eine Wahrscheinlichkeitsverteilung festlegen. Schwierig wird es aber, zumindest sobald wir keine ganz simplen Modelle mehr haben, bei dem Teil P(X), der Wahrscheinlichkeit der Daten. Hierzu müsste man die Wahrscheinlichkeit der Daten über alle überhaupt möglichen Parameterwerte integrieren und das ist oft nicht möglich. Zum Glück kann dieses Problem aber mit einer Simulationstechnik, die in den 80er und 90er Jahren entwickelt wurde, umgangen werden:

6.2.1 MCMC - Markov Chain Monte Carlo

MCMC-Verfahren können Wahrscheinlichkeitsverteilung annähern (in diesem Fall P(H|X)), die man nicht analystisch berechnen kann. Die intuitivste und ziemlich geniale Erklärung, die ich dazu gefunden habe ist die von Michael Franke auf seiner github Seite https://michael-franke.github.io/intro-data-analysis/Ch-03-03-estimation-algorithms.html Ich habe sie hier in Teilen übersetzt:

Image by brgfx on Freepik

Wie die Äpfel an die Bäume kommen
Jedes Jahr im Frühling schickt Mutter Natur die Kinder los, um Äpfel in den Apfelbäumen zu verteilen. Für jeden Baum soll die Anzahl der Äpfel proportional zur Anzahl der Blätter sein: Riese Ralf mit doppelt so vielen Blättern wie Dünner Daniel soll also am Ende auch doppelt so viele Äpfel haben. Wenn es also insgesamt \(n_a\) Äpfel gibt, und \(L(t)\) die Anzahl der Blätter von Baum \(t\) ist, soll jeder Baum \(A(t)\) Äpfel bekommen.

\[A(t) = \frac {L(t)}{\sum L(t')}n_a\]

Das Problem ist nur, dass Mutter Natur \(\sum L(t')\) (die Normalisierungskonstante) nicht ausrechnen kann, also nicht weiß, wie viele Blätter alle Bäume zusammengenommen haben.
Die Kinder (Markov-Ketten) können aber zählen. Sie können zwar nicht alle Bäume besuchen und sich die Zahlen auch nicht so lange merken, aber immerhin für je zwei Bäume. Was sie jetzt tun, ist folgendes: sie starten je an einem zufälligen Baum (Parameterwert), hängen schonmal einen Apfel rein, zählen dort die Blätter \(L(t_1)\) und suchen dann nach einem Baum in der Nähe, von dem sie auch die Blätter zählen \(L(t_2)\) (der Zähler unserer Verteilung). Ist die Zahl der Blätter im zweiten Baum höher, gehen sie dort auf jeden Fall hin und hängen einen Apfel hinr ein. Ist sie niedriger gehen sie nur mit der Wahrscheinlichkeit \(L(t_2)/L(t_1)\) dorthin und hängen einen Apfel auf. So laufen sie weiter zwischen den Bäumen hin und her und am Ende sind die Äpfel ungefähr richtig verteilt. Je mehr Äpfel es sind, desto besser das Ergebnis. Die Besuchshäufigkeit der Kinder hat sich also der gewünschten Verteilung (die Mutter Natur nicht ausrechnen konnte) angenähert!

MCMC Verfahren machen das gleiche: eine MCMC Kette zieht zufällig Werte für die Modell-Parameter und berechnet damit Ergebnis1 = Likelihood der Daten * Prior. Dann zieht sie zufällig Werte, die um die ersten Werte herum liegen, und berechnet dafür wiederum Likelihood der Daten * Prior = Ergebnis2. Wenn Ergebnis2 höher ist als Ergebnis1, springt sie dorthin und zieht von dort aus neue Parameter-Werte. Wenn das Ergebnis2 niedriger ist, springt sie nur mit der Wahrscheinlichkeit Ergebnis2/Ergebnis1 dorthin. Jetzt werden wieder zufällig Werte gezogen usw.
In der folgenden Abbildung sind die erfolgreichen Sprünge durch blaue Pfeile symbolisiert, die abgelehnten Sprünge durch grüne Pfeile.

Visualisierung einer MCMC-Simulation. Sie müssen sich vorstellen, dass Sie die lachsfarbene Parameterlandschaft, die die Werte für Likelihood * Prior darstellt, nicht sehen können, sondern nur die Werte der einzelnen Punkte berechnen können. Weitere Erklärungen im Text. Quelle: https://www.turing.ac.uk/research/research-projects/adaptive-multilevel-mcmc-sampling

Wenn dieser Prozess lang genug fortgesetzt wird, nähert sich die Verteilung der blauen Punkte der a-posteriori Wahrscheinlichkeitsverteilung der Parameter an, die wir haben wollen: Wir haben a-priori Wissen und die Information, die in den neu gesammelten Daten steckt, erfolgreich zum a-posteriori Wissen kombiniert!

6.3 Ein (etwas konstruiertes) Beispiel

Wir müssen uns entscheiden, welche Hühnerrasse wir in größerem Maßstab auf unserer Streuobstwiese halten möchten. Wir haben selbst schon 3 Kraienköppe und 6 Niederrheiner. Die Niederreihner sind uns ein bisschen lieber, weil sie weniger aggressiv sind, aber die Kraienköppe scheinen eine etwas höhere Legeleistung zu haben. Ist es deshalb wirklich sinnvoll, eher Kraienköppe zu wählen?

Schritt 1: Definition der A-priori-Wahrscheinlichkeitsverteilung

Wir erkundigen uns über die Legeleistung der beiden Rassen und erfahren, dass Niederrheiner zwischen 190 und 210 Eier pro Jahr legen, Kraienköppe zwischen 200 und 220. Entsprechend formulieren wir unsere a-priori Wahrscheinlichkeitsverteilungen:

library(brms, warn.conflicts=F, quietly=T)

## Loading 'brms' package (version 2.23.0). Useful instructions
## can be found by typing help('brms'). A more detailed introduction
## to the package is available through vignette('brms_overview').

# mit der Funktion set_prior definieren wir 2 Normalverteilungen, eine mit Mittelwert 200 
# und Standardabweichung 1,5, und eine mit Mittelwert 210 und Standardabweichung 1,5. 
# Unter 'coef' legen wir fest, zu welchen Parameter im Modell, das wir später 
# fromulieren, diese Priors gehören.
priors <- c(set_prior("normal(200, 5)", coef = "rasseN"), set_prior("normal(210, 5)", coef = "rasseK"))

#So können wir sie plotten, um ein Gefühl dafür zu bekommen
curve(dnorm(x, 200, 5), from = 170, to = 230, xlab="Legeleistung", ylab="Dichte")

2. Erhebung von neuen Daten

Jetzt geben wir unsere eigenen Beobachtungen der Legeleistung der 3 Kraienköppe und 6 Niederrheinern aus dem vorigen Jahr an:

rasse <- c("K", "K", "K", "N", "N", "N", "N", "N", "N")
ll <- c(225, 222, 221, 219, 219, 216, 221, 218, 217)
daten <- data.frame(rasse=rasse, ll=ll)

3. Kombination der a-priori-Verteilung mit den Daten, um die a-posteriori-Verteilung zu erhalten

Um die a-posteriori Verteilung zu schätzen nutzen wir die Funktion ‘brm’ des Packets ‘brms’. Wie wir es schon von anderen Auswertungen kennen, formulieren wir zuerst unser Modell, nämlich das die Legeleistung ll von der Rasse abhängt. Die -1 in der Modell-Formel bewirkt, dass das Modell die Mittelwerte für die Legeleistung schätzt und nicht nur den Mittelwert für die erste Rasse und den Unterschied dazu. Unter ‘data’ geben wir unsere selbst erhobenen Daten ein und unter ‘prior’ die oben definierten Priors. Das argument ‘silent’ bestimmt, wie viele Informationen auf der Konsole ausgegeben werden, wenn die MCMC Ketten loslaufen. Diese Simulationen laufen je nach Modell-Komplexizität einige Sekunden oder Minuten.

legeleistung <- brm(ll ~ rasse -1, data = daten, prior = priors, silent = 2)

4. Interpretation

Der obige Code hat die MCMC Simulation durchgeführt und wir können uns jetzt die a-posteriori Verteilungen für die Legeleistungen der Hühnerrassen anschauen:

# So kann man sich die Zusammenfassung anschauen
summary(legeleistung)

##  Family: gaussian 
##   Links: mu = identity 
## Formula: ll ~ rasse - 1 
##    Data: daten (Number of observations: 9) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Regression Coefficients:
##        Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## rasseK   221.74      1.50   218.21   224.29 1.00     1991     1238
## rasseN   217.60      1.16   214.98   219.42 1.00     2180     1756
## 
## Further Distributional Parameters:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     2.35      0.92     1.30     4.63 1.00     1313     1654
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Die Zusammenfassung zeigt uns zuerst nochmal unsere Eingaben und ein paar Informationen zu den Markov-Ketten. Am interessantesten für uns sind natürlich die Schätzer für rasseK und rasseN und ihre Kredibilitätsintervalle. Wie Sie sehen, überschneiden sich die Konfidenzintervalle: der untere Wert von rasseK (l-95%) ist mit etwa 218 kleiner als der obere Wert (u-95%) von etwa 219 von rasseN, die legeschwächer ist.) Sigma ist die geschätzte Standardabweichung der angenommenen Normalverteilung und interessiert uns hier weniger.

plot(legeleistung)

In der Abbildung sind die Ergebnisse noch leichter zu interpretieren. Links sehen Sie die a-posteriori Verteilungen für beide Rassen (beachten Sie, dass die x-Achsen der beiden oberen Verteilungen nicht ausgerichtet sind). Rechts daneben die Dynamik der Markov-Ketten. Aus den a-posteriori Verteilungen ließe sich jetzt auch die genaue Wahrscheinlichkeit berechnen, dass die Legeleistung wirklich unterschiedlich ist. Da wir aber schon gesehen haben, dass die Kredibilitätsintervalle sich überschneiden, bleiben wir bei unserer Vorliebe für die Niederrheiner, auch wenn es nicht ausgeschlossen ist, dass sie ein paar Eier weniger legen.

6.4 Frequentistisch versus Bayesisch

Die sinnvollere Interpretierbarkeit der Ergebnisse von bayesischen Analysen wird am besten durch folgene Abbildung dargestellt:

Quelle: Korner-Nievergelt und Hüppop (2016). Fünf mögliche Resultate der Schätzung eines Effektes, zum Beispiel der Unterschied im Ertrag bei anderer Düngung. Die Punkte geben die geschätzten Differenzen, die vertikalen Balken die 95 % Unsicherheitsintervalle (Vertrauensintervall bzw. Kredibilitätsintervall) an. Die Resultate des zugehörigen Nullhypothesentest sind in der ersten Zeile aufgeführt. In der zweiten Zeile finden sich die posterior Wahrscheinlichkeit für die Hypothese, dass der Effekt „wirtschaftlich relevant“ ist. Die Hintergrundfarbe unterscheidet schematisch „wirtschaftlich relevante“ (orange) von „wirtschaftlich nicht relevanten“ (hellblau) Effektgrößen.

Bei einem unkritischen frequentistischen Ansatz würde das Ergebnis lauten: “Der Effekt bei Gruppe/Behandlung A und E ist signifikant, also Düngen wir jetzt mit dem getesteten Mittel. Bei den anderen Gruppen gibt es anscheinend keinen Einfluss.”
Natürlich gibt es auch in B einen wichtigen Effekt und auch wenn der Effekt in E statistisch signifikant ist, ist die Effektgröße auf den Ertrag so klein, dass es von wirtschaftlichem Standpunkt aus wahrscheinlich keinen Sinn macht, zusätzlich zu düngen. Die Ergebnisse der bayesischen Analyse (also die Wahrscheinlichkeiten für die Hypothese, dass es einen Effekt gibt), spiegeln diese Interpretationen viel direkter wider.

Auch in diesem Kapitel gibt es eine Quiz-Fragen.